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ABSTRACT 


An idealized steady state model of a stream of energetic electrons 
neutralized by a reverse current in the pre-flare solar plasma is 
developed. These calculations indicate that, in some cases, a signifi- 
cant fraction of the beam energy may be dissipated by the reverse current. 
Joule heating by the reverse current is a more effective mechanism for 
heating the plasma than collisional losses from the energetic electrons 
because the Ohmic losses are caused by thermal electrons in the reverse 
current which have much shorter mean free paths than the energetic 
electrons . 

Analysis of the steady state model indicates that it can not 
adequately describe the interaction of the beam with the solar plasma 
because the atmosphere is rapidly heated. If the time scale lor this 
heating is short enough, the density of the atmosphere can be taken 
constant in time. The charge separation required to drive the reverse 
current is expected to respond to changes on a time scale very short 
compared to the time for the ambient plasma temperature to change signi- 
ficantly, so it is a reasonable approximation to use the steady state 
results lor the electric field. With these simplifications, the heating 
due to reverse currents is calculated for two injected energetic electrai 
fluxes. For the smaller injected flux, the temperature of the coronal 
plasma is raised by about a factor of two. The larger flux causes the 
reverse current drift velocity to exceed the critical velocity for the 
onset of ion-cyclotron turbulence, producing anomalous resistivity and 
an order of magnitude increase in the temperature. The heating is so 
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rapid that the lack of ionization equilibrium may produce a soft x-ray 
and EUV pulse from the corona. 
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1. INTRODUCTION 


This dissertation examines the consequences of reverse ciu’rents that 
may be expected to develop in the solar atmosphere in response to the 
imposition of a directed stream of energetic (non-thermal) electrons. 

The phenomena which indicate the presence of streams of electrons mani- 
fest themselves primarily in the "flash phase" of solar flares (Svestka 

19750* flares exhibit a "flash phase" (Svestka 197?^ Sweet 19o9; 

Sturrock and Coppi 19b6) and the existence of directed streams of non- 
thermal electrons is not universally accepted ( Svestka 197!?, Brown 197^; 
Bi'own and Melrose 1977)* A short historical review is presented (cf, l.l) 
as an attempt to place the phenomena in pex’spective. Obsei’vations that 
indicate the pi’esence of energetic electi'ons in the solar atmosphere ai-e 
reviewed and the introduction concludes with a short summary of our 
present theoretical undei’standing of the flare process. 

In Chapter 2 the objections to unneutralized electron beams and 
previous wox*k on reverse curx'ents ai-e summarized and a steady state model 
of a stream of energetic electrons neuti’alized by a revei'se current is 

developed. In Chapter 3 the model is modified to include time dependence 

for a restricted case. The results of Chapters 2 and 3 are summarized 
in Chapter ^1- and possible extensions of the present work are suggested. 
Details of the numerical calculations of Chapters 2 and 3 are discussed 
in the appendices. 

1.1 Historical Overview 

The svm is the closest star to the earth and the only star which we 
can presently observe in great detail. Aside from the intrinsic interest 
of solar phenomena^ we can hope that by understandiTig solar phenomena we 
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will gain insight into what is likely to happen on other stars like the 
sun. The sun is a normal G type main sequence star, but by virtue of its 
position it is the brightest object in the sky. The importance of the 
sun to life on earth cannot be overstated. In the introduction to his 
book. The Sun, C. A. Young (l902a) emphasizes this point, 

"it is true from the highest point of view the sun is 
only one of a multitude - a single star among millions - 
thousands of which, most likely, exceed him in brightness, 
magnitude and power. He is only a private in the host of heaven. 

"But he alone, among the countless myriads, is near 
enough to affect terrestrial affairs in any sensible degree; 
and his influence upon them is such that it is hard to find 
the word to name it; it is more than mere control and dominance. 

He does not, like the moon, simply modify and determine certain 
more or less important activities upon the surface of the 
earth, but he is almost absolutely, in a material sense, the 
prime mover of the whole. To him we can trace directly nearly 
all the energy involved in all phenomena, mechanical, chemical 
or vital. Cut off his rays for even a single month, and the 
earth would die; all life upon its surface would cease." 

The great preponderance of the energy flux from the sun is, to the 

best of our knowledge, very nearly constant (Smith and Gottlieb 197^). 

It is only in those portions of the electromagnetic spectrum where the 

solar output is small (radio XIP/, X-ray), in individual spectral lines 

(e.g. Ca H and k), and in particle emission (the solar wind, energetic 

electrons and nuclei), that the sun’s output varies significantly due to 

solar activity. 

The most obvious manifestations of solar activity are sunspots. 
Sunspots have been observed telescopically since 16H, shortly after the 
invention of the telescope, and with the unaided eye on infrequent occa- 
sions since ancient times (Bray and Loughhead 1964). It is not clear 
which of four men, Galileo Galilei, Johann Goldsmid, Thomas Harriot or 
Christopher Scheiner, actually made the first telescopic observation of 
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svinspots (Bray and Loughhead 196 ^) . That another manifestation of solar 
activity, faculae, were observed at about the same time is demonstrated 
by the title of Christopher Scheiner's ( 1630 ) book, Rosa Ursina Sive Sol 
ex Admirando Facularum and Macularum Fuarum Pheonomeno Varius ( see Eddy 
et al. 1977 ^ Meadows I97O) . lu the fii’st half of the 19th century Schwabe 
( 1844 ) announced the possible existence of the sunspot cycle with a per- 
iod of about 10 yeai's ("von ungefahr 10 Jahren”), Wolf (1852) later 
deduced a more accurate period of 11.1111 ± .038 years or "de sorte que 
neuf periodes equivalent Justement a un siecle". Wolf (iBpS) also deduced 
from eai’lier recoi’ds the years of svinspot minima back to I7OO, but the 
earlier portion of this histoi'ical reconstx'uction has been questioned 
recently (Eddy I976) . 

The first recorded observation of a solar flare occui'i'ed on September 
1 st, 1869 - ^ relatively rare "white light flare", visible against the 
photosphere, was simultaneously observed by Carrington (1859) 

Hodgsoix ( 1869 ) • Id 1808 Janssen (l 869 ) ^nd Lockyer (I869) independently 
discovei’ed that prominences could be seen outside eclipse with a spectro- 
scope with a wide entrance slit. Thereafter various obsei'vers, espe- 
cially Secchi (1877) extensive visual observations of the fonns of 

the chromosphere and prominences using this technique. Flares in 
ixvdividual lines were obsei'ved quite often fi'om the 1870's onward (see 
Young 1871, 1902b, c foi' early examples). The first photographs of flares 
were obtained by Hale (I892) with a spectroheliograph of his own inven- 
tion (Hale 1891). Deslandres ( 1893 ) independently developed a similar 
instrument, and the basic pi'inciple of the spectroheliograph was known to 
Janssen (I869) who actually constructed an instrument similar to the 
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spectrohelloscope (Millochau and Stefanik, I 906 ) for observing promi- 
nences but abandoned it in favor of a widened spectroscope slit. The 
basic principle was independently discovered by Braun, and Lohse attemp- 
ted the construction of a spectroheliograph (Hale I906) . The matter of 
who actually used a "spectroheliograph" first was the subject of some 
debate between Deslandres and Hale (Hale I 9 C 6 , Deslandres, I 905 ) 
this distinction is generally given to Hale. In I908 Hale (I908) made 
the first observation of magnetic fields on the sun, and realized soon 
thereafter that magnetic fields, sunspots and flares were Intimately 
connected (Hale I 929 ) . Because the spectroheliograph took a relatively 
long time to form an image of the whole sun, the systematic investigation 
of flares did not begin until the spectrohelloscope, constructed by Hale 
in 1926 (Hale 1929)^ fully developed (Smith and Smith I 963 ) . The 

development of the polarizing monochromatic filter (Lyot filter) by 
Oilman in 1938 ( Oilman 1938), independently of Lyot's original proposal 
(Loyt 1933 ^ Evans 19^9) ^ allowed photographs of the entire solar disk in 
one spectral line to be made rapidly. This type of filter is still 
widely used in flare patrol telescopes and solar observatories. 

Jansky (1933) made the first observation of radio emission from an 
extra-terrestrial source. It was not until 19^2 that Hey (1946) dis- 
covered meter wavelength radiation from the sun. At about the same time 
Southworth (1945) discovered centimeter wavelength radiation from the 
sun. Reber (1944) made the first published report of radio emission 
from the sun; the earlier work was not published due to its association 
with the war effort. Appleton (1945) published evidence for radio 
emission from the sun in the 7~30 meter wavelength band. Appleton's 
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results were based on amateur i-adio opei-ators' reports (dating from 
1936 ) of "hiss" heard only during the dayliglit hours and frequently 
before sudden fade outs. Appleton and Hey (lQiR-\) noted that some radio 
bursts were associated with flares, Covington (1948) first reported 
microwave bursts from the sun near the maximvim of solar cycle I 8 . 

Burnight (194-9) I’eported the first observation of X-ray emission 
from the sun. Burnight ’s observation was made using photographic film 
with aluminum and beryllium filters flown in a captured VS rocket. 

Peterson and W’inklei' (19'”'8) made the first observation of a flare asso- 
ciated impulsive X-ray burst using a balloon borne proportional counter. 

1.2 Review of Obsei'vations 

The presence of energetic electrons in the solar atmosphere is 
inferred from impulsive hard X-ray bursts, impulsive microwave bursts 
and observations of energetic electi'ons by satellites in earth ox-bit. 
Impulsive micx-owave bursts are rapid enhancements of i-adio flux at 
fi-eqxiencies gi-eatei- than ~ 1 GHk. These impxilsive enhancements occur 
simxiltaneously with impulsive X-i-ay and EUV bursts and often show very 
similax' time sti'ucture, even down to the fine details of the time profiles 
(Petei-son and Winklei- lSX>9x Kundxi iQi-'l, Andei'son and Winklei- 19o2, Kane 
and Domxelly 1971; deFeitex’ 1977; Svestka 197‘-) • impulsive micro- 

wave bursts are genex-ally attributed to gyro-synchrotron x-adiation from 
electx'ons with enex'gies gx-eater than 100 keV (Holt and Ramaty 19t>9; 
Svestka 1979) • The gradual post-burst increases can be intex-preted as 
thermal bremsstrahlmxg from the flare-associated soft X-ray plasma and 
are usxially accompanied by x-adio eralsslcn at lowex- fx-eqxiencies (Svestka 
1977 ) • The appax-ent discx'epancy between the nxanbex- of electx-ons x-eqxiired 



to produce the impulsive X-ray bursts and the number of electrons required 
to produce the impulsive microwave emission (Peterson and Winkler 1959) 

can be resolved if the details of the microwave emission in the solar 

atmosphere are considered (Holt and Ramaty 1969^ Takakura 1972). The 
generation and propagation of microwaves in the solar atmosphere during 
a solar flare are complicated processes involving the magnetic field 
configuration, ambient plasma density and temperature and density and 
energy spectrum of the non-thermal electrons (Holt and Ramaty 1969^ 

Kruger 1972, Takakura 1972^ Svestka I 975 ) . Therefore, it is difficult to 
unambiguously infer the number and energy spectrum of the non-thermal 
electrons from the observed microwave emission. 

In some flares, non-thermal electrons escape into the interplane- 
tary medium and are observed by satellites in earth orbit ( Svestka 1975 ^ 
Lin 197 ^) • Since the electrons apparently propagate primarily along 
magnetic field lines in the interplanetary medium, electrons are observed 
primarily from flares in the western half of the visible hemisphere of 
the sun or from flares beuind the west limb of the sun (Svestka 1975^ 

Lin 197 ^)- Lin (197^0 concludes that there are two distinct types of 
non-relativistic electron' bursts (e < 5 OO keV) observed at 1 AU, "pure 
electron events", that is those not accompanied by energetic (> 10 mev) 
protons, and "mixed events" during which both energetic electrons and 
protons are observed. The energy spectra of the "pure electron" events 
can be well fitted between 5 keV and 100 keV by a power-law in energy, 
dN/dE “ E , with Y ~ 2.5~5*5 Lut exhibit a rapid steepening at ener- 
gies above 100 keV (Lin 197^)* On the other hand the typical spectra of 
energetic electrons for "mixed events" extend smoothly in a power-law 
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out to and beyond 10 inov and tend to be somewhat harder (Y ~ 

(bln 19'T^I-) . When impulsive X-ray bursts are associated with electrons 

n A 

observed at 1 aU, 10‘"“10 more electrons are required to produce the 
impulsive hard X-ray bursts than escape to the intoi’planetaxy medium 
(bin and Hudson 19'?^^ l*lu 19‘f^l-) . 

It is now genoraiiy believed that the mechanism for the production 
of impulsive hard X-ray bui'sts is bromsstrahluug from electrons scattering 
on protons and heavier ions (Kane Bvostka 19'f!b Brown I9‘f‘0 • 

Smaller impulsive events generally consist of ono or a few spikes with 
comparable e-folding rise and fall times of ->^10 s (Kane Kano and 

Anderson I9'f0j Cranncll ot ai. 19’f'f ) . bargor events, with total durations 
of minutes o,r tens of minutes, usually have a complex spiky time structure 
(Svestka l97’b Hoyng al . l9*(’o) ’ Kx’ost and Dennis and I'rost 

(l9‘f^0 have also reported an apparently distinct non-impulslve nou- 
thormal hard X-ray component in some lai'ger events, after the impulsive 
phase of the flare and possibly associated with a second phase of particle 
acceleration, In this work, we restrict ovir attention to the impulsive 
hard X-ray bursts, and assume that both the later "second phase" hard 
X-rays and the "gradual components" in the low energy channels (<. '>0 keV) 
of some instruments are distinct piionomena. 

The spectral information on impulsive hard X-ray bursts is limited, 
but most events can be reasonably fitted to a decreasing power-law in 
photon energy between 10-20 keV and 100“1‘'0 ke? (Kano lO'fii-, Brown lO'f'b 
Hojmg 19‘rf) • power-law index is typically between ;*!,*' and (Kane 

Svostka lO'}';') although some bursts have very soft spectra and 
power-law indices as largo as have boon repox'tod (Peterson et al. . 
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Most events show a softening of the spectrum at higher energies (Kane 
1974 ^ Svestka 1975) • This bend or "knee" in the power-law spectrum 
usually occurs between 60 and 100 keV (Brown 1975^ Svestka 1975) l>u.t in 
some events can occui’ as high as 5 OO keV (Brown 1975) • Since the high 
energy cut-off of many instruments is below 5 OO keV [e.g, OSO -7 
(Peterson et al. I 974 ), OSO-5 (Frost ^ 1970) OGO -5 (Kane and 

Anderson I 970 )] such a break in the spectrum may be present in many 
events for which no break is reported. It is obvious that the power- 
law must flatten at low energies, otherwise the total X-ray flux would 
diverge. However, the determination of the low energy cut-off is diffi- 
cult because the X-ray emission at low energies (< 10 keV) is dominated 
by the gradual quasi -thermal component in most events (Brown 1975^ 

Svestka 1975) • 

Although the Interpretation of the X-ray spectrum as bremsstrahlung 
from a non-thermal (i.e. non-Maxwellian in energy) distribution of 
electrons is widely accepted, some workers advocate a thermal inter- 
pretation for many impulsive X-ray bursts (for example Chubb 1970^ 

Elcan 1976 ^ Crannell et al. 1977) some events seem to fit an 
exponential rather than a single power-law spectrum (Elcan 197^^ 

Crannell e^ al. 1977)- However, the spectral data are poor, particularly 
at higher energies (primarily due to counting statistics), and it is not 
clear that an exponential spectrum is to be preferred over two power- 
laws or some other form for the spectra. Brown (1974) has demonstrated 
that any observed hard X-ray spectrum can be produced by a thermal 
plasma with a suitable temperature distribution in the source. Brown 
( 1975 ) also pointed out that the emitted X-ray spectrum is rather 



insensitivo to the sovu'ce electi'oiv energy spectrum npci concludes ttat a 
power-law electron spectriim is not strongly mandated by the presently 
available data. 

There are theoretical objections to multi -thermal models of 
impulsive hard X-ray bursts (Kahler and models that produce power 

law spectra by superposing different exponoivtial spectra seem somewhat 
contrived to this author despite assertions to the contrar 5 >^ by some 
workers (e.g. hrown lyf ') < In more recent "thei'inal" models of impul- 
sive hard X-ray bursts (e.g. Smith and Lilliequist the electron 

distribution is not expected to be Maxwolliaix. Neither the theoretical 
objections (Kahler, limited Qbservnvtional snpport for thermal 

electron distributions (Elcan 19 'fo^ Cx'aivnoll et 1977) are relevant 
to this type of model. 

‘fhore is some support, from observations of impulsive EUV bursts, 
for the view that the impulsive hard X-ra,y bursts are produced by non- 
thermal, energetic electrons streaming from the corona to the chromo- 
sphere, Impulsive EW bursts have been observed directly by satellites 
(for example Kelley and Kense l9T9j Hall I 97 I) ai\d indirectly from the 
ionosphorie effects they produce (Doivnelly W'rl) , Those bvirsts show 
close tlmo coincidence with the impnlsivo X-ray and microwave bursts and 
the time profiles closely resemble the X-ray and microwave bursts 
(Noyes 19' (’*!•> Donnelly l';57h‘, Kane and Donnelly l^Tfl, Kane 197*1-) • I'he 
energy radiated in the 10-1030 A band is ^ 10^''-10' times the energy 
radiated in the associated impxulsive hard X-.vny burst (Donnellj*- 19'f*!-, 

Kano and Donnelly iD'fl) , This ratio of exxergy radiated ix\ the EHV and 
hard X-ray bands correspoixds qualitatively with the expected ratio of 



Coulomb colliiDional losses to bremsstrahlung emission from a thick- 
target hard X-ray source (Koch and Motz 19‘39^ Petrosian 1973^ Donnelly 
Brown 1979 ) . There are indications that the P,W radiation origi- 
nates low in the chromosphere. The density of the EUV emitting region 

1? -3 

can be estimated to be ^ 10 cm (Donnelly 197^^ Kane and Donnelly 
1971 ), corresponding to the solar chromosphere, and the EUV bursts 
exhibit (statistical) limb darkening which would be expected if the 
radiation originates in the chromosphere (Kane and Donnelly I 97 I) ■ 
Although the observations presently available do not exclude other’ 
interpretations, the preponderance of evidence seems to favor bremsstrnh- 
lung from a non-thermal distribution of eirergetic electrons as the 
source of impulsive hard X-ray bursts (Svestka I 971 ) . 

If the emergent X-ray spectrum were known exactly, the spectrum of 
the energetic electrons that produce the radiation, averaged over the 
source, could in principle, be recovered (Bi’own 197!') • The two extreme 
approximations that are usually considered are "thick-target" and "thin- 
target" (Brown 197!9 Svestka 197!^j Hudson 197^0* thin-target 

approximation the electrons lose a negligible amount of energy in the 
hard X-ray source (Brown 197!b Svestka 197!'', Hudson 197^) • this 

approximatioii, the mean electron source spectrum [i,e, the instantaiieous 
average of the electron energy specti’um over the emitting volume 
weighted by the background density, see Brown (1979)] is just the spect- 
rum of accelerated electrons. In the thick-target approximation, this 
is not the case. 

In the thick-target approximation, the electrons lose all their 
energy (primarily by Coulomb collisions) in the source regioii. Since 
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the mean free path of the more enei'getic electrons is longer^ the energy 
spectrum of the electrons averaged over the emitting volume is hardei'. 

If we assume the density is »mifoi*m in the source^ that the electrons 
are all streaming downward and that the accelerated electron enei’gy 
spectrum is a fairly steep power-law, then the approximate difference 
between the inferred mean electron source spectxnim in the thin and 
thick-target approximations can be estimated. In this case, since the 
mean free path of an electron against Coulomb collisions is approxi- 
mately pi'oportional to , the effective source volume lor electrons 
of energy E in the accelerated spectrum is also approximately propor- 
tional to 1^ . Since the injected specti'um is very steep (by assumption), 
once an electron has lost an appreciable fraction of its enex'gy, it no 
longex’ contributes significantly to the emex'gent X-ray flux. Therefore, 
to produce the same power -law index of emex'gent X-x-ays the index of the 
injected electx-on beam must be ~ g greater (a softer injected spectx'um) 
in the thick-target case than in the thin-target approximation. The 
pX’eceeding simple nnalj^sis neglects beaming effects in the case the 
exiex'getic electron velocity distx'ibution is anisotropic (Petrosian 19?3^ 
Brown 1978 ) ^nd the exact behavior of the Coulomb cross section. However, 
the conclusion is found to be qxialitatively cox'rect in thick-tax-get 
models of impulsive hard X-ray bxirsts for X-ray enei'gies below ~ 100 keV 
even when a more detailed analysis is performed (Brown 197?^ Petrosiaxx 
1973 ^ Hudson 1972 ^ Brown 1971)- The more detailed calculations Indicate 
that, depending on the assumptions and model characteristics, thick- 
tax'get models reqxiire Injected electx'on power-law indices I. 9-2 greater 
than thin-target models fox' the same emex'gent X-x'ay spectra. 


11 



Some early models of impulsive X-ray bursts considered impulsive 


injection of the energetic electrons and the subsequent decay of the 
impulsively injected electrons in the source region (e.g. Takakura and 
Kai 1966). In its simplest fonn this model does not agree with observa- 
tions since it predicts a systematic hardening of the burst spectra 
during the decay of the hard X-ray burst (Brown 1975; Petrosian 1973) 
for the same reason that the source averaged energetic electron spectx’um 
is harder in the thick-target models, i.e. the low energy electrons lose 
their energy moi’e rapidly than the high energy electrons. Browxi (1972) 
has introdxiced a modification of the usual coronal impulsive hard X-x*ay 
soux’ce that removes this pax’ticulax* objection, but it requix-es the 
assximption of an average soxu’ce density than is energy dependent 
(n E , o: > 3 / 2 ) . Brown (1972) motivates this assumptioxx by invoking 
an enex’gy dependent pitch angle distribution for the accelerated elec- 
tx'ons, but the simplicity of the original "cox'onal tx’ap" model is lost. 
This type of model caxi explain the obsex’vatioxx of impulsive X-ray bux’sts 
fx’om "behind-the-limb" flax-es since pox'tions of the X-x'ay source ax’e 
high ixx the cox’oixa. Howevex’, since behind-the-limb flares also prodxice 
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higii energy X-rays, this model requix*es densities ^ 10 cm high 

, 9 \ 

(^10 cm) in the solar atmosphere. If this were the case, impulsive 
X-ray bxirsts from behiixd -the -limb flares coxild be explained by thick- 
target models as well. Since the product of the instantaneous number' of 
energetic electrons in the source and background density determines the 
emergent X-ray flux. Brown's (1972) model x’equires a much larger number 
of energetic (> gO kev) electrons than equivalent thick -target models. 
Additionally, since most of the enex’gy resides in the low energy electx’ons 
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which encounter only low densities 10 ) , these electrons cannot be 
Invoked to account for the impulsive EUV bursts which are emitted from 

12 “3 

regions where the density is ^ 10"^ cm (Donnelly 197^) which are 

observed simultaneously with impulsive X-ray bursts (this is also true 
of more recent "thermal" models, e.g, Smith and Lilliequist I 978 ) . 

Aside from some difficulty in accounting for impulsive X-ray bursts 
from behind-the-limb flares, thick-target models for the hard X-ray 
bursts are at least not excluded by present observations. Since they 
have the advantage of also providing the energy required lor the impul- 
sive EUV bursts (Donnelly 197^)^ it seems I'easonable to accept the thick 
target approximation lor the production of the hard X-ray bursts. In 
this case, since the time for the electrons to lose all their energy is 
short compared to the time scale of the impulsive X-ray burst (Brown 
1975 ^ Petrosian 1973), variations in the X-ray flux and spectrum are 
attributed to changes in the (unspecified, c.f. 1.3) acceleration process 
In the thick-target model the energy flux of the electron stream 
required to produce a specified X-ray flux at 1 AU depends on (a) the 
anisotropy of the electron velocity distribvitioia, (b) the power-law 
index of the X-ray flux and (c) the lowest enei’gy to which the power-law 
in energy is assumed to extend for the enei’getic electrons (Brown 1975, 
Petrosian 1973) • Neglecting possible beaming of the bremsstrahlung 
radiation (petroSian 1973, Brown 1972) and backscatter from atmosphere 
(banger and Petx’osian 1977)> obtain an order of magnitude esti- 

mate for the flux of non-thermal electrons at the sun for an observed 

flux of X-rays at 1 AU. If the flux of X-rays at some energy at 

-2 -1 -1 

eai'th is J (photons Cm s keV ), then the total X-ray photon flux 
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is 10 ' y photons keV s . Since the observed power-law spectra 

are typically fairly steep (Brown 1975^ Svestka 1975)# most of the X-rays 
at Eq are produced by electrons with only slightly higher energies. The 
total efficiency (the ratio of bremsstrahlung losses to Coulomb collisional 
losses) is approximately 10~^E for a thick-target hydrogen plasma 
(Koch and Motz 1959# Petrosian 1973). Therefore, the total non-thermal 
electron flux in the source above E^ must be 

1.3 Review of Flare Theories 

In the preceeding sections we have discussed some of the observed 
properties of solar flares as they relate to the inferred presence of 
non-thermal electrons in the solar atmosphere during a flare. We have 
not dealt with most of the diverse phenomena associated with solar 
flares. Svestka (1975) lists thirty-seven "basic properties of flares". 
Wheii all the subtopics are counted, Svestka' s list contains more than 
eighty observational aspects of flares. With such a large number of 
properties to be considered, it is not surprising that a wide variety 
of flare theories and models have been proposed. Since reviews of flare 
theories exist in the literature (Svestka 1975# Sweet I 969 )# the selec- 
tion of theoretical ideas discussed here is only representative and not 
exhaustive. This discussion of flare theories is included only to show 
how the production of non-thermal electron streams fits in the present 
theoretical picture of solar flares and therefore no particular model 
will be treated in detail. 

It is now widely believed that the energy released in solar flai’es 
is stored in the magnetic fields in the upper solar atmosphere (Rust 
1977 # Svestka 1975# Sweet I 969 ) * enei’gy which is available for 
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release Is the excess energy of the non-potential magnetic field con- 
figuration above the energy in the potential (current-free) field (Rust 
19TT> Svestka 1975> Sweet I 969 ) • Because the magnetic field energy 
density is generally believed to be greater than the thermal energy 
density of the plasma in the upper solar atmosphere, the non-potential 
field configurations must be nearly force-free (Gold and Hoyle I 96 O, 

Sturrock 197^) • Although many non-potential field configurations have 
been proposed, these configurations can be divided into two broad cate- 
gories depending on the distribution of currents in the solar atmosphere 
(Svestka 1975^ Sturtock 197^)* One possibility is a force-free configu- 
ration in the form of twisted flux tubes (Gold and Hoyle I 96 O, Alfven 
and Carlqvist 1967 ^ Spicer 1977) sheared field lines (Tanaka and 
Nakagawa 1973) . In this case the curreiats are distributed over a large 

;■ ■ I 

volume in the atinosphere. The other possibility is that the field is 
largely current-free With the current concentrated in current sheets . 

(Sweet 1959^ Syrovatsky I966, Sturrock I968, Priest and Heyvaerts 197 ^)- 
A large number of flare models have been developed under the assumption 
that current sheets develop in the solar atmosphere as a result of 
motions in the photosphere or the emergence of new flux (svestka 1975 ^ 

Sweet I9S9) • Barnes and Sturrock (I972) have studied the development of 
non-potential force-free fields due to photospheric motions and found 
that the stored energy in the force-free configuration can exceed that 
of a configuration with a current sheet. They concluded that one possible 
sequence of events that would produce a current sheet in the solar atmo- 
sphere was the conversion of a more energetic force-free configuration 
to a configuration With a sheet. Priest and Heyvaerts (197^-) examined 



the pi’oduction of a current sheet when new flux emerges Into a pre- - 

existing magnetic field, configuration. 

The earliest electi*omagnetic models of flares invoked the production 
of non-thenual electrons and realized the importance of electi’ic fields 

at "neutral points" in the magnetic field (Giovanelli 19 'h5, 1948, 

1 

Hoyle 1948 ) . Dungey (19!58) pointed out that, when reconnection of 
Magnetic field lines occurs, a DC electi'ic field will be developed in 

I ■' ■ ;■ 

the reconnection region which could lead to acceleration of chai’ged 
particles. In "cui*a*ent interruption" models (Alfven and Carlqvist ,19 oT) 
ele<ibrons are accelerated by the DC electric field that develops when 
the "inductive circuit" is opened. In models in which reconnection 
occurs in a current sheet (sturi-cck 1988 , Fx’iedman and Hamberger 1969 # 
Coppi and Fi’iedland 3 , 971 ), some acceleration by a DC electric field at 
the neutral point may occur, but the bulk of the acceleration is usually 
attributed to stochastic acceleration of electrons by high frequency 

. . I ' : . 

elecjtric fields that develop during the recomxection process due to 
plasma instabilities (sturrock 1974,* Smith 1974). It has proved diffi- 
cult to develop a self consistent theoretical model of the rapid 
acceleration of the number of electrons required to produce the observed 
X-ray flux (Smith 1977^,* Brown and Melrose 1977)- present, the 
mechanism by which electrons are accelerated in the impulsive phase of 
solar flares is not well understood theoretically (Svestka 197;'0 • 

HO|Wever, simple considerations indicate that if the energy stored in the 

] 

magnetic field is released in the low density corona, particles can be 

. I 

expected to acquire enex'gies of lO-lOO key ( Sturrock 1974)^ Furthexmiore, 
the ingredients of many possible accelex'ation mechanisms (DC electric 
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fields, plasma turbulence) are natural by-products of most processes 
which release the energy stored in the magnetic field. Therefore, 
since; there is observational evidence for the acceleration of electrons 
in' the impulsive phase of flares, we will assume that this acceleration 
does occur even though the exact mechanism has yet to be elucidated. 



2. STEADY STATE MODEL OF BE.^I AND REVERSE CURRENT 


2.1 r Objections to Unneutralized Beams 

The simplest thick-target model for the production of impulsive 
Xi-ray bursts is that considered by Petrosian (1973) • In this model a 
beami of energetic electrons is assumed to propagate from the corona to 
the chromosphere. All the electrons are assumed to have their velocities 
in the same direction until they lose all their energy, approximating a 

source in which the energetic electrons stream down a nearly vertical 

I 

magnetic field line with small pitch angles into an atmosphere with a 
small density scale height (Petrosian 1973). Several authors (Brown 
1976 ^ Brown and Melrose 1977^ Colgate e^ al . 1977> Hoyng 1977^ Hoyng ^ 
al. 1976 ) have pointed out difficulties if this electron stream is not 
neutralized by a reverse current. 

! Brown (1976) pointed out that the number of electrons required to 
stream from the corona to the denser portions of the solar atmosphere 
during some impulsive hard X-ray bursts was quite large. Indeed in some 

events as large as 10 ^^ (Hoyng et al. I976), or all the electrons in the 

. 2-3 - 

solar atmosphere above the level where the electron density is 10 cm 

(fiBrown I 976 ). Another objection to the existence of an unneutralized 
beain is that the magnetic energy that would be stored in this beam is 
many orders of magnitude larger than the total flare energy (Colgate et 
al. 1977 ). If N is the total number of electrons streaming downward 

Over the duration t(s) of the impulsive phase, the magnitude (emu) of 

I 

the current may be estimated from 

I =- ec'^NT 



( 2 . 1 ) 



I If the transverse and longitudinal dimensions of the stream are of 
order R (cm), an estimate of the strength B (gauss) of the magnetic 
field produced by the stream is given by \ 

i B =* 2IR“^ ; (2.2) 

and the total energy U (ergs) of this magnetic field may be estimated 
from 

I 

1 ^ ? 1 9 

u - ^R-^B - — I^R , (2.3) 

which becomes 

U e^c“^l^T“^R ' (2.4) 

K^e and Anderson ( 1970 ) estimate the total energy involved in a 

pq 2 

typical small flare to be ~ 10 ^ ergs, the time scale to be 10 s , and 

8.5 

the characteristic length scale to be 10 ’ cm and infer from the X-ray 

35 

data that the total number of energetic electrons is 10 . For these 

13 2 5 

values the above formulae lead to estimates of I /«=< 10 ' , B 10 , 

34 

and U (^10 . For a large event the total flare energy could be 

32 9 5 3 

^ 10~^ ergs, the length scale 10"^ the characteristic time 10 

iud the total number of energetic electrons «=< 10^^ (Hoyng et al. I 976 ). 

16.2 7 4l 

In this case I -v 10 * , B pa 10 and U 10 . Clearly a model which j 

involves an unneutralized beam leads to unacceptably high values of the 
magnetic field and magnetic energy associated with the beam. 

Problems associated with the propagation of high current beams of 
charged particles not neutralized by a reverse current have been considered 
in other contexts. Alfveh (1939) examined the limitations on the propa- 
gation of electrostatically neutralized high current beams of relativistic 

19 , 



charged particles, motivated by an apparent sidereal day variation in 
the cosmic ray flux (Alfven 1938, Compton and Getting 1935) ^ which 
later proved to be spurious (Dorman 19T^) . Consider a cylindrically 
symmetric, mono-energetic, uniform beam of charged particles moving 
through a background of opposite charge so distributed that the charge 
density (esu) is everywhere zero. If the beam is infinite in extent 
along the symmetry axis and has a radius of R , then the magnetic field 
as a function of distance from the axis for r ^ R is 

B(r) = = 2Jtjr , (2.5) 

where l(r) is the current inside r and j is the current density 
(assumed uniform). The gyro radius of a charged particle in a magnetic 
field is 



( 2 . 6 ) 


where p is the particle momentum and q is the particle charge. 

Consider a test particle of the same charge and mass as the beam particles 
moving in the magnetic field of the beam. Suppose the test particle is 
initially at the outer edge of the beam (i’=R) and has the same momentum 
as the beam particles. We denote by R^ the beam radius for which the 
gyroradius of this particle in the average field it sees in its trajectory 
is equal to the beam radius. For a beam of this radius (R^) ) the 
particle will cross the axis of symmetry with its momentum perpendicular 
to that of the beam particles. If the radius of the beam is increased, 
the particle will cross the s 3 Tnmetry axis with the component of its 
momentum opposite in sign from that of the beam particles and its average 
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velocity over the trajectory will also be negative. Clearly increasing 
the beam radius beyond will not increase the current. If we esti- 
mate the average magnetic field as 1/2 the field at the edge of the 

beam, we find that there is a maximum current which can be carried by a 
beam which satisfies our original assumptions ; 


Therefore 



(2.7) 



o X/2 

here p=v/c and Y=(l-P ) . I is called the "Alfven current limit" 

J\ 

or the "Alfven-Lawson current limit" and for electrons we find 


- 1700 PY , 


(2.9) 


in agreement with Alfven' s (1939) more rigorous derivation. This restric- 
tion is much more stringent than the objections to the stored magnetic 
energy. Tor an electron energy of 100 keV, the currents estimated for 
the hypothetical small and large events are 10^^ and 10^^ 

respectively. The value of the current limit derived by Alfven depends 
on all the original assumptions being satisfied. Arbitrarily large 
currents can in principle be propagated by relaxing the assumption of 
exact electrostatic neutralization (Lawson 1957^ I 958 , 1959)> 
assumption that the beam is mono-energetic (Bennett 193^)^ the assumption 
that the current density (particle flux) is uniform (Hammer and Rostoker 
1970 ) adding a very strong magnetic field along the symmetry axis 
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(Hammer and Rostoker I97O) . However, none of these mechanisms seem 
particularly likely to be applicable in solar flare impulsive hard 
X-ray bursts, although some are relevant to particular laboratory experi- 
ments. The simplest resolution to these objections is the existence of 
a reverse current (cf. 2.2). 

2.2 Previous Work on Reverse Currents 

It is well known that a plasma tends to preserve charge neutrality. 

A process which tends to give an excess positive or negative charge in 
some region will lead to electric fields which act upon the plasma. 
Movement of electrons in response to this electric field will then 
restore charge neutrality. One expects that analogous process will also 
tend to maintain current neutrality. If an electron beam is suddenly 
introduced into a plasma, a sudden change occurs in the magnetic field 
structure which will develop induced electric fields opposing the primary 
current. 

Although interest in beams of relativistic electrons is not recent 
(see for example Bennett 193^, Alfven 1939)^ theoretical and experimental 
work on high current relativistic electron beams was stimulated by the 
development of devices capable of producing relativistic electron beams 
with currents on the order of or greater than the Alfven-Lawson current 
limit (see for example Graybill and Nablo I966, Roberts and Bennett I968, 
Yonas and Spence I969) • Roberts and Bennett (I968) injected a beam of 
3.5 nisv electrons (p=. 992 ^ Y=T* 85 ) with a beam current of 3000 emu 
(l *23l.) into a linear pinch with 10 '^cm . They found that 

the beam current was nearly completely neutralized by a reverse current 
in the ambient plasma and that the change in the total current (measured) 
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was a very small fraction of the beam current. Similar results have 
been obtained with other experimental apparatus (prono e^ al. 1975^ 

Ekdahl et al. IS'jh, Goldenbaum e^ al . 19 T^> Klok et al, 197 ^^ Miller and 
Kuswa 1973, Levine et al, I97I) when the ambient plasma density was 
sufficiently high. 

Several theoretical models of energetic electron beams neutralized 
or partially neutralized by reverse currents in the ambient plasma have 
been developed (lor example Cox and Bennett 197 ^^ Hammer and Rostoker 
1970, Lee and Sudan 197^^ Lovelace and Sudan 1971^ Chu and Rostoker 1973) 
Since these theoretical . treatments are primarily concerned with the high 
current energetic electron beams that are typically produced in labora- 
tory studies and not in the electron beams thought to be responsible for 
impulsive hard X-ray bursts, some of the results are not relevant to the 
solar flare case (cf. 2.3). The models cited treat cylindrically symmet- 
ric mono-energetic beams of the type considered by Alfven (cl, 2.1) with 
the possible addition of a uniform magnetic field along the symmetry 
axis. When the beam current is small compared to , then the induced 
reverse current flows primarily outside the beam cylinder (r > r) while 
for I » the reverse current is confined to r ^ R and the current 
neutralization is local in the sense that the ambient electrons drift 
with the velocity 



( 2 . 10 ) 


where V, is the reverse current drift velocity, V, is the velocity of 
d b 

the beam electrons and n, and n are the beam and plasma electron 

be 

number densities (Cox and Bennett I97O) . Depending upon the sharpness 
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of the leading edge of the beam, large amplitude coherent plasma oscilla- 
tions may be generated by the passage of the beam head (Hammer and 
Rostoker 1970, Cox and Bennett 1970^ lee and Sudan 197^^ Ohu and Rostoker 
1973 ) . The amplitude of these plasma oscillations is t) ^ , where 

e /m^) j and T is the rise 

time of the beam (Lee and Sudan 1971)* These oscillations decay with a 

scale length of , where is a phenomenological momentum 

relaxation time for the plasma electrons. If the lateral dimension (r) 

of the beam is large compared to the electromagnetic skin depth 

(X =c/o) ) , after the plasma oscillations decay the net current will be 
£ p 

~ X /R times the beam current. The current of the beam will be neutral- 
E 

ized for a length of V, T (R/X„) ^ The theoretical models for mono- 

b ee' E 

energetic beams are not appropriate for the streams of energetic electrons 
that are responsible for impulsive X-ray bursts. We argue in Section 2.3 
that the beams in solar flares will be current neutralized in steady 
state. 

2.3 Steady State Model 

We now examine a simple model for an impulsive X-ray burst. We 
consider a vertical flux tube extending from the corona to the chromo- 
sphere and assume that electrons are accelerated at the top of the flux 
tube by the development of stochastic electric fields (sturrock I966, 

Hall and Sturrock 1967^ Newman 1973) or by some other mechanism (cf. 
Section 1.3)» The injection of these electrons down the field toward the 
chromosphere then leads to the development of a reverse current both by 
the mechanisms considered for mono-energetic beams in laboratory plasmas 
(cox and Bennett 197^^ Hammer and Rostoker 1970^ Chu and Rostoker 1973) 
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and due to an electrostatic field due to charge imbalances. The strong 
tendency of a plasma to remain charge neutral implies that, if a current 
is generated in the plasma that would systematically violate dp/dt=0 
on time scales much greater than a plasma period (i.e. a non-MHD current) 
then this current will generate a neutralizing secondary reverse current. 

In contrast to the mono-energetic beams typical of laboratory experi 
ments, the streams of energetic electrons that produce impulsive X-ray 
bursts probably have smooth distributions in energy. This is inferred 
from observations (cf. 1.2) and theoretical considerations indicate it 
is likely that the number of electrons does not increase with energy 
(Brown and Melrose 1977^ Smith 1975) • consider below an energetic 
electron stream with a distribution of this type, that has electrons of 
all energies present. The low energy electrons are constantly merging 
with the background plasma and can build up charge Imbalances. In the 
case of a mono-energetic beam considered by other workers (for example 
Cox and Bennett 197^; Chu and Rostoker 1973) > charge imbalance would only 
build up at the ends of the plasma device since the energetic electrons 
do not interact with the plasma significantly except through the reverse 
current. Charge built up at the ends of an experimental plasma column 
would either be conducted away by external return paths or be shielded 
from the bulk of the plasma within a few Debye lengths of the ends and 
not drive reverse currents in most of the volume of the plasma column. 

Lovelace and Sudan (I 97 I) pointed out that the microscopic process 
involved in heating the plasma with reverse currents are equivalent to 
heating with currents induced by external fields. However, the reverse 
currents avoid the skin effect limitations of currents induced by 
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external fields. Similarly, since charge can be supplied by the beam in 
the case of solar flares, charge imbalances can build up within the plasma 
and drive reverse currents. Although these charge imbalances arise 
throughout the plasma, we can estimate the time (t^) required to accumu- 
late sufficient charge separation from the time required to accumulate 
enough charge per unit area on a parallel plate capacitor to produce an 
electric field sufficient to drive the required reverse current. This 
required chai'ge separation is related to the current density by 


jT) = E il-it S = Uxj c T , (2.11) 

' unn c 

where £ is a sui’face charge density, T| is the resistivity, and 
is the unneutralized portion of the beam current density. Then the time 
to accumulate the requii’ed charge is 



( 2 . 12 ) 


The ratio of unneutralized cui'rent density to the beam current density is 
\^/R (cf. 2.2) so that 


■^c = 


IV 

knc^ 


1/2 


T ^ 10 ^•^^(loV)(n/10^) (R/10^) . 


(2.13) 


(2.14) 


This assumes that the resistivity is the usual Spitzer value. If the 

resistivity is "anomalous" the effective collision frequency can be of 

order the electron plasma frequency (;o ) . Actually this is an upper 

P 

limits for the Buneman instahility the effective collision frequency is 



'V .lai (Buneman 195^) • resistivity is proportional to the collision 

P 

frequency so we can write 


Ti ln^*'75„V2 

~ , (2.15) 


%pitzer 10^'^^nT 


SO that T becomes 
c 


T 

C 




(2.16) 


We see that the time to accumulate charge imbalances sufficient to drive 
a neutralizing reverse current is short compared to time scales of 
interest. 

If the resistivity is written 


m c ^ 

^ e 1 

^ ~ 2 T ’ 

n e ee 

e 


(2. IT) 


then T becomes 
c 


T = — (ai T ) 
c c p ee 


( 2 . 18 ) 


and the ratio of the charge accumulation time to the time the current 
remains neutralized (“^^^) 'th® mechanisms considered for a mono- 

energetic beam (cf. 2.2) becomes 


R / \ —1 

— c ^'^p ee^ _ / 

T R R / N “ ^^P ^ee^ ; 

n - - (o) t ) 

c X ^ p ee' 


(2.19) 


= 10“^-^^(a, (1 o9/R)(1oV)^'^^ , (2.20) 

n 
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so that the charge accumulation time is much shorter that the time the 
current would remain neutralized if no charge imbalance arose in the 
plasma. For time scales and length scales of interest in solar impulsive 
X-ray bursts, the reverse currents will be caused primarily by charge 
separation (Hoyng and Melrose 1977) • Also, since beams of interest for 
solar impulsive X-ray bursts are not expected to have sharp fronts, the 
plasma oscillations excited by passage of the "beam head" will be of 
extremely small amplitude and consequently of no great significance 
(Melrose 197^) • Therefore, we are justified in considering a steady 
state in which the beam current is exactly balanced by a reverse current 
in the background plasma. For the present (cf. Chapter 3)^ we assume 
that the background plasma can be adequately described by a Maxwellian 
velocity distribution and use transport coefficients based on this 
assumption (Sptizer I 962 ) . 

We are interested in the case in which, the primary electron stream 
is composed of high-energy electrons with consequently long mean free 
paths in the tenuous solar corona. However, we shall find that the 
electric field that develops to drive the reverse current also decelerates 
the electron stream (cf. Lovelace and Sudan I 97 I) • But when the electron 
energy becomes comparable with the thermal energy, the mean free path 
will be sufficiently short that the primary electrons will merge with 
the background plasma. As a simple representation of this process, we 
ignore collisions in discussing tJae primary beam but we assume that an 
electron of the primary beam is absorbed into the background plasma 
when it is decelerated to zero energy. This approximation is justified, 
if the temperature of the ambient plasma is sufficiently low. 
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as a further simplification, we consider a flux tube of uniform 
cross section, we may use the following simple one-dimensional form of 
the Vlasov equation: 


§f 

a 


e d(f> Sf 
m ds iv 


( 2 . 21 ) 


Where s measures length along the tube, v is velocity (along the tube), 
f(s, v) is the velocity distribution function of the primary electron 
stream, and 4’ is the electrostatic potential. 

At the top of the flux tube ( s=0) , the primary electron stream is 
moving with positive velocity and electrons that are decelerated to zero 
velocity are assumed to be removed from the beam. Hence we may without 
ambiguity, express f in terms of 'F , which is defined by 

. ( 2 . 22 ) 

The initial distribution function may therefore be expressed as 

f(0,v) = f('F) . (2.23) 

With this initial condition, we find that the solution of the Vlasov 
equation (2.21) is 


f(s, v) = F(f-4>) . 


(2.24) 


The current density in the primary electron stream is given by 


^s = - c 


J Hs, 


'■)v dv , 


(2.25) 


which may be expressed as 
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(2.26) 


me J 
0 


F(Y-4>)d'K . 


Since will prove to be negative in the region of interest, it is 
convenient to write 


G = -4> , 

so that Equation (2.26) may be reexpressed as 


(2.27) 


00 

---/ 

*0 


F(x)dx 


(2.28) 


We have seen that the beam current will be nearly completely neutralized 
by currents in the background plasma, so we may write 


^p + js = ° ^ 


(2.29) 


where is the secondary current induced in the background plasma. We 

here assume that the density and temperature are such that j may be 

P 


represented by Ohm's law. 


jp = = 1] 


-1 ^ 
ds 


(2.30) 


It is convenient to introduce a new independent variable | to 
replace s by the relationship 

= T| ds (2.31) 

Then, on substituting Equations (2.28) and (2.30) into Equation (2.29) 
and differentiating with respect to ^ , we obtain 
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(2.32) 


— F(9) ^ = 0 . 
^^2 me d5 


It is convenient to solve this equation foi' 5 terms of 

? = X(G) , 

rather than vice versa. Equation (2.32) becomes 

^ - V(Q) 

9.[aQj ,q2 - ^ ^ 


which may be integrated once to give 


me 


9 

=/*F(Q' 

ide) J J 


)dG' , 


if we assume that 9=0 (cj)=0) and X = 0 (^=0) at s = 0 . We 
from Equations (2.28), (2.29) ^iid (2-30) that 


e 

me 


F(0')dG' 


F=0 0 


Hence Equation (2.35) becomes 


dX _ me 

dG ~ 2 
e 


E(9')d0' 


Le 


It is now convenient to inti'oduce a specific foimi for F(f) 


f(y) = K(¥q+'F)‘ 


© . 

(2.33) 

(2.34) 


(2.35) 

find 

(2.36) 


(2.37) 

(2.38) 
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This is a power-law distribution at high energy which flattens at low 

energy, with the "knee" characterized by . 

We introduce the symbol h(T, s) for the flux of electrons 
-2 -1 

(cm s ) of energy exceeding ef at the position s : 


eo 

^s) = f(¥' + Q(s))dY' 


(2.39) 


If the initial flux is written as H (f) , we find that 


HnOf) ' 


so that the total particle flux is given by 


' ivfjir ■ (2.41) 

With the form of Equation (2.38) for F('if) , Equation (2.3T) 
integrates to give 

X(9) = ^ ('!' + G)^ ~ '‘'ol . (2.h2) 

We easily obtain from Equation (2.42) an expression for the (negative) 
electric potential Q in terms of the resistivity weighted distance 
measure ^ : 

2 

e(C) = (''o yT Sr s) - »o • (2-‘*3) 

Hence from Equation (2.39)^ find that 
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■ 0 

me / 


-(Y-l) 




On noting that the electric current carried by the stream is related 
to h(Y, s) by 


J<,(s) = - I h(0,s) , 


(2.45) 


we see that 


jgCs) 


2 / 2 
e^K / Y ^ Y e^K 

(Y-l)mc I 0 Y-l me 


-{(Y-l)/Y] 


( 2 . 


In order to specify the current, particle flux, and electric field 
as functions of s , we must adopt a specific form for i|(s) . A con- 
venient approximation to the density and temperatui’e structure of the 
solar atmosphere, which is expressible in analytic form, is px'ovided by 
the constant heat flux model. If we now assume that s measures distance 
vertically downward from the corona, and that n = n^ and T = at 

s = 0 , this model (Adams and Sturrock 1975) yields the following 
expressions : 


t(s) = , 


(2.47) 


n(s) = nQ[TQ/T(s)] exp {- 


/ 7/2 bFs^5/7 ^/21 

(Tq -bFs) - ^ Jj, (2.11-8) 


- 1.21 6.58 -2 -1 

where a 10 ' , b 10 , and F (ergs cm s ) is the downward 

heat flux. 

The resistivity, in modified Gaussian units, may be derived from the 
expression given by Spitzer (I 962 ) : 
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(2.U9) 


. T ^ 3.64 

where g 10 

related to s by 


T1 = 

Hence we find from Equation (2.31) that 5 is 



(2.50) 


Our model is then completely specified by the choice of the coronal 
temperature, the coronal density, the coronal heat flux, Y , the energy 
corresponding to , and the injected energetic electron flux. For the 
coronal parameters, we adopt values typical of the corona above an 
active region (Noyes I 97 I) : 

T 3 X 10^ K , 

9 -3 

n 10-^ cm , 

6 -2 -1 
F ^ 5 X 10 erg cm s 

We choose to correspond to 25 keV ; and we choose Y =2,5 ^ The 

fraction of the beam energy deposited and the total energy deposited by 

6 4 

Joule heating between T = 3 X 10 K and T = 3 x 10 K as a function 
of the energetic electron flux are displayed in Figure 2.1. For a flare 

1.9 *5 2 

area of 10 cm ^ the energetic electron flux inferred from a large 

It “2 “1 

impulsive X-ray burst corresponds to 3^ 10 ' cm s (Hoyng et al. I 976 ). 

Figure 2.2 illustrates the energy deposition rate due to Joule heating 

as a function of temperature of the atmosphere for this injected energetic 

electron flux. The ordinate of Figure 2.2 is the time required to raise 

7 

the ambient plasma temperature by 10 K, if the plasma were heated at the 
steady state rate. As we will see (cf. Chapter 3)^ the heating rate 
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NUMBER FLUX (crn“2sec“') 


ENERGY DEPOSITED (erg cm'^sec'*) 



Figure 2.1 The fraction of beam energy deposited (solid curve) and the 
total energy deposited (broken curve) by Joule heating as a function of 
the energetic particle number flux. 



10° 1.5x10° 6x10°. 3x10° 1.5x10° 6x10^ 3: 

T(K) 

Figure 2.2 Joule heating rate as a function of ambient plasma temperature 

17 -2 -1 

for an energetic particle number flux of 10 cm s . The heating rate 

7 

is displayed as the time T required to heat the plasma by 10 K . 



decreases as the temperature of the plasma increases, so the ordinate of 
Figure 2.2 is only representative of the heating rate immediately after 
the beam is turned on. 

\Ve can check the assumption that Coulomb collisions are not impoi'- 

tant for the energetic electrons in the beam. Since the reverse current 

losses are proportional to the resistivity for a constant current density, 

-3/2 

these losses are proportional to T , The Coulomb collisional losses 

for energetic electrons of the same kinetic energy are proportional to 

n . Therefore, we expect the ratio of reverse current losses to co] - 

3/2 -1 

lisional losses to be proportional to (nT ) . Reverse current 

losses will be more important than Coulomb collisional losses for an 
energetic electron in the beam if 

a = io"‘^'\v/v^)^(v^/v^) > 1 , (2.51) 

where V is the velocity of the energetic electron, V is the electron 
thermal velocity, and (cf. Equation 2.10) is the reverse current 

drift velocity (Hoyng et al. 197^) • expected, for the same kinetic 

3/2 -1 

energy and current density, o: is proportional to (nT ) since 

1/2 -1 

T and V^j n . if we define the injected energetic electron 
flux as the flux of electrons with kinetic energies greater than eY^ 
then from Equation (2.4o) we find 

= • ^ 2 - 52 ) 

(Y-l)m2^ ^ ° 

Since the reverse current drift velocity is related to the current 
density by 
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V = £ 

d en 


(2.53) 


we may write the drift velocity in terms of the injected energy flux as 


V 


d 



. -[(Y-1)/Yl 


( 2 . 5 ^) 


For the adopted values Y=2.5 and eY^=S5 heV, we find that the ratio 
of reverse current losses to Coulomb collisional losses (n) for a 
beam electron with kinetic energy 25 keV in the adopted constant heat 
flux model atmosphere is 


a 


10 ' 


,2.82 




nT 


3/2 


1+10 


-23.56 


F 


0 


-[ (y-i)/y] 


(2.55) 


In Figure 2.3, the energy at which a=l is plotted as a function of 
temperature for several values of H . We see that for any energetic 
electron flux we have considered, the energy at which Coulomb collisions 
are as impor1:ant as the reverse current losses for the energetic electrons 
is reasonably low, indicating that oui- assumption that Coulomb collisions 
may be neglected is an adequate approximation. 


38 



KINETIC ENERGY (keV) 



3x10® 1.5 


6 3x10® 1.5 

T(K) 


Figure 2.3 The kinetic energy for which the reverse current losses are equal to 
Coulomb collisional losses for a beam electron as a function of temperature in 
the constant heat flux atmosphere, for five different energetic particle number 
fluxes. From top to bottom the curves correspond to energetic particle number 


fluxes (cm^s of 


and 10 



3. REVERSE CURRENT HEATING 


3.1 Generalization of Steady Model to Time Dependent Case 

As we have indicated, for the energetic electron fluxes required to 
account for the observed X-ray flux by thick-target bremsstrahlung, the 
ambient plasma is rapidly heated by the reverse current. The rate at 
which the background plasma is heated by the reverse current depends on 

the beam current density and the ambient plasma density and temperature. 

1/2 

If the ratio of V to the electron thermal velocity (v = (2kT /m ) ) 

d ^9 ^ ^ ^ 

is large enough, the background plasma may be unstable to the growth of 
electrostatic plasma turbulence which can dramatically enhance the plasma 
resistivity and therefore, the reverse current heating rate. For example, 
the reverse current will be unstable against the excitation of ion- 
acoustic or electrostatic ion-cyclotron turbulence unless for T^ and 
T^ the electron and ion temperatures, respectively, (Kindel and Kennel 

T^ — .1 T^ (ion-acoustic turbulence) 

T^ ^ .3 T^ ( ion-cyclotron -turbulence) 

Te T^ (ion -cyclotron turbulence) 

T^ 3 ^1 (ion-cyclotron turbulence) 

T^ =* 10 T^ (ion-acoustic turbulence) 

Reverse current heating is a self quenching process. If the reverse 
current is stable against the growth of electrostatic turbulence, then, 
as the plasma is heated by the reverse current, the resistivity decreases 
and the reverse current losses are reduced. If the reverse current is 
unstable to the growth of electrostatic turbulence, the plasma will be 


19T1) 


d t, e 


2-5 for 
.9 for 
^ .3 for 

.1 for 
.05 for 


1+0 



heated until the instability criterion is no longer satisfied. The 
heating of the plasma will also cause a pressure imbalance. The time 
T (s) for the plasma to .respond to changes of pressure by bulk motions 
can be estimated from 


T-L/V^^. , (3.1) 

Where L is a characteristic length and . is the ion thermal velo- 

t, 1 

7 

city. Even for a temperature as large as 10 K, this time is long 

2 9 • T 

(lO s) compared with the heating time for a length scale of 10"^’ cm, 
so that the plasma density will not change appreciably during the 
heating. Since we expect reverse currents to be established locally on 
time scales on the order of a plasma period 10 ^ s) which is much 
shortei- than the time scale fox’ heating of the plasma 10 s), it 

should be a reasonable approximation to use the results of Chapter 2 
for the instantaneous velocity distribution of the energetic electrons 
as a function of distance from the injection point. 

We have calculated the heating due to reverse currents for two 
Injected energetic electron fluxes (h ) , The heating rate was taken 

ill 

to be jxist that which results from the Olimic losses suffered by the 
reverse current and is given by 


3t _ c 
at " 3n^K 



(3.2) 


The electron and ion temperatures were assumed to be equal. We shall 
say more about this assumption later. At each time step the current at 
each spatial grid point was calculated using Equation (2.46). The time 
step was regulated so that the largest change in temperature at any grid 
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point was 1^ in one time step. Since we have not found an analytic 
solution to the time dependent problem considered here, the constant 
heat flux model of the atmosphere was abandoned in favor of a more 
accurate numerical model which is discussed in Appendix B, The spatial 
grid spacing was chosen so that for the initial temperature profile 
(t=0) the temperature change between spatial grid points was less than 
1'^. The atmosphere was assumed static; that is, the number density (n) 
was held constant in time. The details of the numerical methods used 
are discussed in Appendix A. 

We have used the same Y and as in Chapter 2. The results 

Iv 

for an injected energetic electron flux of l,4l4 x displayed 

in Figure 3*1 while similar curves for an injected energetic electron 
flux of 5*656 X 10 are displayed in Figure 3*2. Figure 3*3 depicts 
the density structure of the model atmosphere. The abscissa, I , of 
the figures is integrated number density from the injection point, 
defined by 

s 

n(s')ds‘' , ( 3 . 3 ) 

0 

where n is the total number density (sum of neutral hydrogen and 

proton density) . Because we have used a nuirterical model rather than 

the simple analytic constant heat flux model, we Were not free to choose 

the density at the injection point (see Appendix b) . The initial density 

in the adopted model is approximately twice the density in the constant 

heat flux model used in Chapter 2. Since the reverse current heating 

2 

rate is proportional to j and inversely proportional to density, the 




Figure 3.1 Temperature (t) as a function of Integrated 
number density (l) from the injection point, for an 
energetic electron number flux of 1.4l4 X 10^"^ cm~^ s~^. 
The temperature is displayed before the beam is injected 
(t = 0 s) and for two times after the beam is injected 
( t = 1 s and t = 4 s ) . 



I (cm"^) 


Figure 3.2 Temperature (t) as a function of integrated 
number density (l) from the injection point, for an 
energetic electron number flux of 5.656 X 10 ^"^ cm~^ s~^. 
The temperature is displayed before the beam is injected 
(t = 0.25 s and t = 1 s). 





I (cm^) 


Figure 3.3 Number density of neutral and ionized hydrogen 
(n) as a function of integrated number density (l) from 
the injection point. The model serves only to represent 
the gross overall structure of the solar atmosphere above 
an active region (see Appendix b) . 




lower energetic electron flux corresponds roughly to the initial heating 


rate shown in Figure 2.2. 

Figure 3.1 shows the temperature as a function of I for two times. 

Is and 4s after the injection of the beam. The energetic electron flux 
used in the calculation of rhe results displayed in Figure 3.2 is four 
times that used for Figure 3.1. Figure 3-2 displays the temperature 
after . 25 s and Is corresponding to the same total energy input as the 
curves lor Is and 4s in Figux’e 3-1. Thermal conductivity was neglected 
in these calculations, but computer runs with thermal conductivity 
included indicated that thermal conductivity did not have significant 
effects for the short time scales 4s) involved in here (see Appendix A) . 

3.2 Anomalous Resistivity and Reverse Current Heating Rate 

The electrical resistivity used irj the calculations depended on the 
reverse current drift velocity as indicated below; 

( ^ 13 . 

I d t, 1 

11 = '! , (3.>1) 

where V . = (2KT./m.) for T. the ion temperature, m. is the 

1 X X X X 

proton mass, is the resistivity due to Coulomb collisions derived 

by Spitzer (I 962 ), and is an anomalous resistivity due to the 

'A 

presence of electrostatic ion-cyclotron turbulence calculated by lonson 
( 1976 ) . For the pui'poses of calculating the value of the anomalous 
resistivity we have adopted B=d-00 gauss, a reasonable value lor the 
pre-flare corona. Since we are considei'ing a flux tube of constant 
ci’oss section, the field is the same for all values of s , the distance 


46 


from the injection point. For the smaller energetic electron flux, the 
reverse current drift velocity did not exceed the critical velocity for 
the onset of electrostatic ion-cyclotron turbulence. In this case the 
temperature of the tenuous coronal plasma was raised by a factor of 2, 
but most of the beam energy was deposited in the dense portion of the 
model atmosphere. 

The larger energetic electron flux, however, caused the reverse 
current drift velocity to exceed the critical velocity for the onset of 
electrostatic ion-cyclotron turbulence in the low density portion of the 
atmosphere, resulting in an anomalous resistivity and an order of magni- 
tude increase in the temperature in these regions in a relatively short 
time. Since Coulomb collisions were neglected in this calculation, the 
heating of the denser portion of the atmosphere is not calculated accu- 
rately after thc; first few tenths of a second (also see Appendix b) . If 
collisions were taken into account for the primary electrons in the 
beam, the heating of the denser regions below the corona would be more 
localized and higher temperatures would be reached. However, these 
results indicate that an energetic electron beam may significantly heat 
the low density coronal plasma much more rapidly than would be calculated 
from considering only the effects of Coulomb collisions on the beam 
electrons . 

The time for electron and ion temperatures to equilibrate by 
Coulomb collisions assuming only one species rather than both species 
are heated as we have assumed may be estimated (Spitzer I 962 ) from 

1 / m \3/2 

t . 12.6 n" T + -^ T. 

ei \ e m. 1 / 

X 


s 


(3.5) 


For the temperatures, densities, and the time scales considered here. 
Coulomb collisions alone will not establish equal electron and ion 
temperatures . We have taken the electron and ion temperatures to be 
equal for computational convenience; however, and must, therefore, 
address the question of whether one species is preferentially heated. 

For the case depicted in Figure 3.1^ for which the resistivity is 
just classical Spitzer resistivity, only the electrons are heated at 
first. According to Equation (3-5)^ th® ions are not likely to be 
heated significantly in turn by energy exchange with the electrons. The 
heat capacity of the plasma is therefore reduced by a factor of 2, and 
the times given in Figure 3-1 should simply be reduced by a factor of 2. 

The situation in which plasma turbulence develops, as for the case 
depicted in Figure 3.2, is considerably more complicated. As we have 
indicated, the critical drift velocity for the onset of electrostatic 
ion acoustic or ion-cyclotron turbulence depends on the ratio of the 
electron and ion temperatures. Just what happens when this drift velo- 
city is exceeded is not well understood, however. 

The anomalous resistivity which we have assumed to result from the 
presence of electrostatic ion-cyclotron turbulence was calculated by 
lonson ( 1976 ) under the assumption that the turbulence saturates by ion 
resonance broadening (Bum and Dupree 19T0) that the electron and ion 

temperatures were equal. Palmadesso et al. (197^)> on the other hand, 
made the first of these assumptions and calculated heatiiig rates of 
electrons and ions. They found that the ions are heated much more 
rapidly than the electrons, and Papadopoulos (197?) subsequently 
concluded that the instability turns off when the ion heating has 
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proceeded to the point at which the instability criterion is no longer 
satisfied. If only the ions are heated, the situation will differ from 
that depicted in Figure 3-2. The temperature plotted should be inter- 
preted as the ion temperature (note that this affects the calculation of 
the expected excitation and ionization rates) and the times given reduced 
by a factor of 2 for the same reason those in Figure 3.1 should be 
reduced if only the electrons are heated. 

It has also been suggested that the ion-cyclotron turbulence satu- 
rates, not by ion resonance broadening, but by the formation of a plateau 
on the electron velocity distribution, instead, in which case no signi- 
ficant anomalous resistivity results (Papadopoulos 197?) • H tliis 
happens, then as in the case without plasma turbulence, only the electx'ons 
are heoted at first, at a rate given approximately by classical resisti- 
vity. In this case, however, larger electx'on beam current densities 
must have been involved to begin with in order for the reverse current 
drift velocity to have exceeded the critical velocity for the onset of 
ion-cyclotron turbulence. Since j is larger than in the case without 
turbulence, the classical heating rate is higher lor this case. If the 
electrons are heated sufficiently in this manner, the critical drift 
velocity for the onset of ion-acoustic turbulence will be exceeded. In 
that case, the electrons will be heated until the criterion for insta- 
bility is no longer satisfied, or until 

y ^ Q (kT /m.) , 

d s ' e i'' ' 

where is the ion sound speed. 


(3.6) 



This latter scenario for rapid electron heating would apply^ for 
instance, to an electron beam strength equal to that assumed in Figure 
3.2. More precisely, assuming the ion-cyclotron turbulence does satu- 
rate by electron plateau formation, a beam of this strength would result 
first in electron heating given approximately by the results in Figure 
3.1 with a time scale reduced by a factor of about 32. After roughly 
•5 s, ion-acoustic turbulence would develop, resulting in rapid heating 
of the electrons to a final tempei'ature which may be estimated from 

T m.V^/k =“ 10^^ K . (3.7) 

In short, the exact behavioi’ of the ratio of the electi’on and ion 
temperatures is not well undei’stood and cannot be determined without a 
much more detailed analysis than is appropriate fox' the present work. 

We have assumed that the electron and ion temperatures are about equal 
as a useful and reasonable approximation with which to estimate the 
magnitude of the x’evex’se current heatiiig. As discussed above, however, 
temperatux’e enhancements much lai'gei’ than those depicted in Figures 3-l 
and 3.2 are possible. 
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k. CONCLUSIONS AND SUGGESTIONS FOR FURTHER RESEARCH 


We have examined a simple model for the production of impulsive hard 
X-ray bursts during solar flares. The model involves a beam of energetic 
electrons propagating from the corona to the chromosphere. We have 
found that if this beam is to exist, the current carried by the beam 
electrons must be neutralized by a reverse current in the background 
plasma. The requirement that the reverse current exist has two conse- 
quences that have not been previously recognized in the context of this 
type of simple model of impulsive hard X-ray bursts. The reverse 
current heats the ambient plasma and the electric field that is developed 
to drive the reverse current decelerates the primary electrons. Joule 
heating by the reverse current is a more effective mechanism for heating 
the tenuous coronal plasma than Coulomb collisional losses from the 
enefgetiq electrons, because the ohmic losses are caused by thermal 
electrons in the reverse current which have much shorter mean free paths 
than do the energetic electrons. 

We have found that the time scale for heating the ambient plasma 
by reverse currents can be comparable with the time scales characteristic 
of impulsive X-ray bursts (Hoyng et al. I976) . It is possible that 
thermal bremsstrahlung from the rapidly heated plasma can account for a 
significant portion of the observed impulsive X-ray flux. Hence this 
mechanism can offer an explanation of the fact that some flares first 
produce high-energy X-ray emission near the top of a loop rather than at 
the footpoints of the loop (Brueckner 19T6) • Anothex- important conse- 
quence of this process is that, if thermal emission can account for a 
substantial fraction of the impulsive flux up to 50 keV, then the 
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number of electrons required to produce the nonthermal X-ray flux is 
greatly reduced (Brown 1975 ) • 

The time scale for heating can also be short compared to the 
ionization times of the plasma ions and may therefore produce non- 
equilibrium line-emission strength enhancements of lines present in the 
plasma spectrum just prior to the rapid heating (Shapiro and Knight 
1978) . These non-equilibrium effects are likely to be observable only 
if plasma turbulence develops causing a large enhancement in the plasma 
resistivity (Shapiro and Knight I978) . 

We have made several simplifying assumptions in ordei' to facilitate 
the calculations presented in Chapters 2 and 3 - In a more realistic 
model, some or perhaps all of these restrictions could be relaxed. We 
now briefly discuss how the relaxation of some of these assumptions is 
likely to change the conclusions we have drawn and suggest possible 
extensions of the work we have presented in Chapter 3 . We have assumed 
that all the electrons in the beam are moving in the same direction, or 
equivalently that they have zero pitch angle. The reverse current 
arises to balance the flux of electi’ons in a given direction due to any 
anisotropy in the energetic electron velocity distribution. If the 
energetic electron velocity distribution is neai’ly isotropic, no signi- 
ficant reverse current will arise (see for example Smith and Lilliequist 
1978) • Even if the distribution is strongly anisotropic, but the 
electrons streaming down from the corona to the chromosphere have non- 
zero pitch angles, the Coulomb collisional losses will be enhanced 
relative to reverse cui'rent losses since the collisional losses are 
proportional to the total path length of the energetic electrons in the 
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atmosphere, while the reverse cui’rent losses are proportional to the 
average component of the energetic electron velocity along the field. 

Since the emergent X-ray spectrum is relatively insensitive to the angular 
distribution of the energetic electron velocities (hanger and Petrosian 

1977) » i't is extremely difficult to infer the reverse current drift 
velocity from measurements of the X-ray flux. A more detailed discussion 
of this and other difficulties in inferring the reverse current drift 
velocity from X-ray observations caai be found elsewhere (Hoyng et al. 

1978) . 

We have neglected the effects of Coulomb collisions on the primary 
electrons. As Figure £.3 demonstrates, this is an adequate approxima- 
tion immediately after the flux of energetic electrons is initiated; 
however, Coulomb collisions become relatively more important as the 
plasma is hegted since the reverse current losses are reduced, Until a 
significant increase in the density of the coronal plasma is effected 
by the evaporation of material from the chromosphere. Coulomb collisions 
are unlikely to be important in the upper portions of the atmosphere. 

In the lower lying dense regions. Coulomb collisions will rapidly domi- 
nate over reverse current losses, and, as we have indicated, affect 
the heating of this portion of the atmosphere. Oiie extension of the 
work presented in Chapter 3 that should provide additional insight into 
the behavior of energetic electrons in the solar atmosphere during flares 
would be to perform a calcvilation similar to that we have presented, but 
include the effects of Coulomb collisions and a distribution of pitch 
angles for the energetic electi’ons. 


We have neglected the dynamics of the background plasma. As we 
have indicated, the rapid heating of the plasma can cause a large pres- 
sure imbalance. Por the results presented in Figures 3.1 a^nd 3.2, the 
pressure is a factor of ^ 20 higher in the high density portion of the 
atmosphere than in the low density regions indicating that evaporation 
of high density material would occur if the dynamics of the ambient 
plasma were accounted for. This would not have a large effect on the 
calculations presented in Chapter 3 because the time scales considered 
are so short. However, on longer time scales mass motions in the atmo- 
sphere could have important effects. Previous work with fluid dynamic 
models of solar flares (for example see Kostyuk and Pikel'ner 1975^ 
Kostyuk 1975; Craig and McClymont 1976) has not included reverse current 
losses. The development of a numerical fluid dynamic model of the 
solar atmosphere heated by a beam of energetic electrons, including 
reverse current losses could provide valuable information about the 
formation of the quasi-thermal soft X-ray plasma that is produced during 
solar flares. 

We have not calculated either the radiation from the heated plasma 
OX’ the bremmstrahlung from the energetic electrons. Since almost all 
the information we now have and are likely to accumulate in the fore- 
seeable fxiture about solar flares comes from the observation of the 
emitted radiation, it would be useful to calculate the emitted radiation 
from any realistic model to ascertain to what degree it resembles the 
solar atmosphere during a flare. 

More realistic models than those we have considered that include 
the effects of Coulomb collisions, the dynamics of the background plasma. 
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a reasonable magnetic field configuration^ radiation and thermal conduc 
tion are necessary to account for the complicated phenomena that are 
observed in solar flares. However^ our study of the reverse current 
and the heating it can cause indicates that reverse currents can play 
an important role, at least in the initial heating of the solar plasma 
during a flare. 


55 



Appendix A 


NUMERICAL AND CCMPUTATIONAL METHODS 

As we have indicated in Chapter 3^ the results for the time depen- 
dent case are calculated by using the steady state results for the 
current as a function of distance from the injection point and calcula- 
ting the change in temperature from a suitably discretized form of 
equation (3*2). In reality, the calculation is done for the more 
general case of partially ionized hydrogen. Since the reverse current 
heating calculation is only accurate in the tenuous high temperature 
portion of the atmosphere, this generalization did not have a substan- 
tial effect on the results of the calculation. However, the manner in 
which the partial ionization is included could in principle be accurate 
in any astrophysical plasma that is sufficiently tenuous that the gas 
is optically thin to its own radiation, photo-excitation and ionization 
are unimportant and collisional de-excitation can be ignored. The 
ionization state of the plasma is then a function of temperature only 
provided non-equilibrium effects can be ignored. The only elements in 
astrophysical plasmas that are sufficiently abundant for their ioniza- 
tion potential to affect the heat capacity of the gas are hydrogen and 
helium. Only hydrogen is included in the present calculation, but since 
the effects of partial ionization on the heat capacity are included via 
a pretabulated interpolation table (discussed below) the effects of 
helium could be included with only minor modification. The modified 
version of equation (3.2) actually solved numerically is 
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where T is defined by 
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T = (1 + X)t + X 
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(A.2) 


where X(t) is the fraction of the hydrogen nuclei that are ionized and 
^lON ionization potential of hydrogen. That is, the thermal 

energy content of the plasma per cubic centimeter is 


"^TH = 2 ’'e 


(A.3) 


The temperature is obtained from T via the interpolation tables men- 
tioned above, and it is obvious that the inclusion of helium only in- 
volves calculating a different interpolation table. In fact we have 
Included only hydrogen and used the expression given by Moore and Fung 

(1972) for x(t) : 


x(t) = + 10 e^[,4288 + i In p + ,4698p"^'^^_ ) 


, (A. 4) 


where p=l5800/t. Then T is implicitly defined as a function of T 

£ 

by Equations (A. 3) and (A.4), 

Spitzer gives the resistivity of a hydrogen plasma as : 


T] = In A , 


(A.5) 


where A is defined by 
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(A. 6 ) 


^ 4.2 X lo^k , 


so that we may write Jj 
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Therefore, T]^ can be written as a sura of a function of T only and a 
function of T only tiraes ln(n) : 


T) = TL(t) + TM(T) ln(n) , 


(A. 8 ) 


where 


tl(t) = 


lo3-31 T 3/2 In T - I In X 


10 ' 
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T S: 4.2 X 10 K , 
(A.9) 


and 


tm(t) = 


lo3.0l ^*3/2 


^q3.24 ^-3/2 


T S 4.2 X lO'^k 


T ^ 4.2 X 10"^k . 


(A. 10) 


The calculation of the current as a function of distance depends 
only on the resistivity weighted distance measure 5 . In Chapter 2 we 
were able to write an analytic expression f oi’ 5 as a function of s , 
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but in the present case the resistivity varies with time. The value of 
I at the ith grid point is approximated by 


5± = 5'i-l + ("Ll * (=1 - ’ 


(A. 11) 


Where superscripts refer to time steps and subscripts refer to spatial 
grid points, and |^ = 0 . So long as the reverse current drift velocity 
is less than the critical velocity for the onset of ion cyclotron turbu- 
lence, the calculation of in this manner is straightforward. How- 

ever, when the background plasma is uijitable to the growth of ion cyclo- 
tron turbulence, the situation is somewhat more complicated. In this 
case the value of 7]'^ depends on the current, and a transcendental 
equation must be solved to find 7]'^ from Equations (2.46), (2.49) and 
the result for anomalous resistivity due to ion cyclotron turbulence 
(lonson 1976 ) 


7]„ = 0.06 (c Q./;» ,)(l - 13V, ./V,) 
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where fl. = (eB/m.c) and :,u = (4rtn e/m ) , we find that is 
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^ 2 n.X^ 
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(A.13) 
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where = cJ'^/enX'? and = 13 V 

d. 1 


t,i 


If we define G(j^) by 


/ Jn j e-K I.V Ye^K 

“(■’i) = ■^1 + r?^yo * ty^ 


me 


5"i-i + 


+ . o6 


' m 


m. 


.1/2 


kne 


(° t - °i- i> 

2 n.X'? 

1 1 


e n.X' 5 ' 
c. 1 i 
1 

cj-^ 


I 

( 


[ (Y-l)/Y] 


(A.ll^) 


then when g(j '?)=0 , satisfies Equation (a, 13 ) • When the current 

calculated neglecting anomalous resistivity corresponds to a drift velo- 
city that is greater than , we take an initial estimate for J'^ ^ 


yl : 


c. 

j'*^ = — - e n. X-? , (A. 15) 

i c 1 X ^ 

and refine this estimate by application of Newton's method to Equation 
(A.li)-) . Examination of Equation (A. 13 ) shows that Newton's method will 
always converge for this initial estimate and the convergence is usually 
reasonably rapid, i.e, usually 6 or fewer iterations are required. 

The functions needed for the calculation [ TL, TM, X and the impli- 
citly defined t(t )] are evaluated by cubic interpolation on pretabula- 
ted tables. The method used is dependent on the architecture of the 
IBM 36O-37G series computers and the internal representation of double 
precision floating point numbers used on these machines. The use of 
pretabulated functions is considerably faster than calls to the FORTRAN 
library routines that would otherwise be necessary. This is particularly 
true in the case of the implicitly defined function t(Tj.) which would 


60 


have to be solved iteratively at each spatial grid point for each time 
step. The internal representation of double precision floating numbers 
on IBM 360-3TO computers is presented diagramatically below, 

hexadecimal exponent (excess 64) 

0 1 'ir 7 8 63 

l4 digit hexadecimal fraction 

sign hit 

FIG. A.l, Internal representation of double 
precision floating point numbers on IBM 360 
and 370 series computers. 

In the ititerpolation procedure^ the first I 6 bits (bits O-lp) are 
extracted, an offset substracted and the result treated as a double 
word displacement from the base address of an interpolation table. The 
remaining 48 bits (16-63) are used to form a floating point fractional 
displacement (frac) from the largest value of the temperature for which 
the function is tabulated which is smaller than the value of the tempera- 
ture for which the value of the function is desired. The value of this 
displacement is such that 0 ^ frac < 1 ; frac is used to calculate 
weights for the four nearest tabulated values of the function in a cubic 
polynomial interpolation. Once the weights fox- the cubic intei’polation 
are calculated only 4 doxible precision floating point multiplies (^ ,6l ps 
each on IBM 370~l68 with high speed multiply; and 3 adds (~ .30 ps each) 
are required to produce an interpolated value from a table. Since the 
weights are to be calculated for TL, TM and X they are also used to 
calculate the critical velocity for the onset of ion-cyclotron tui’bulence. 
This would require, a call to the FORTRAN libi'ary subroutine "DSQRT" which 
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is sufficiently fast that interpolation would be slower for the calcula- 
tion of the square root alone. However^ since the weights must be cal- 
culated for TL, X ^ TM and TE, interpolation is faster than a call to 
"DSQRT" because the weights are effectively "free" for this calculation. 
The semi -logarithmic tabulation scheme allows interpolation from 
T=i|-.096 X 10”^ K 6.710 X 10 K with a maximum relative displacement 
from a value of T for which the function is tabulated of ~ 3?^ with 
only 820 table entries. In fact some of the 820 entries are never used 
due to the nature of the signed magnitude normalized representation of 
floating point numbers on these machines, but the reason for using this 
sort of tabulation scheme is that a reasonably large range can be 
covered with relatively few table entries, and the correct tabulated 
values can be actMssed extremely rapidly. 

Listings of two main programs and several subroutines are provided 
for the sake of completeness. All of the time consuming routines have 
been hand coded in assembly language. Routines that perform initial- 
ization and diagnostic functions as well as the main program are coded 
in FORTRAN. The first main program produces the tables that are required 
for the cubic interpolation. The second main program reads in the inter- 
polation tables, model parameters and starting values. The starting 
values used initially are from the steady state model atmosphere described 
briefly in Appendix B. 

The assembly language programs calculate the current at each spatial 
grid point (CURCAL), calculate the change in temperature at each point 
and determine the time step (tESTP) and write out the arrays at the 
designated intervals (tOUT). In addition the calculation of the current 
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(CURCAL) requires takin.fc the -(y- 1)/Y power of a number, which if done 
with FORTRAN library routines would require taking the natural logarithm 
and exponentiating. Both the library routines "DLOG" and "DEXP" are 
slower than "DSQRT” so an assembly language routine was written that 
calculates the 3/5 power of a number the routine is called by 

CURCAL, The subroutine DTAG is used primarily for monitoring the per- 
formance of the model during program changes and subsequent debugging. 

In production runs it could be replaced with a subroutine that does 
nothing (i.e. returns as soon as it is called) without affecting the 
model calculations; therefore it is not reproduced here. The FORTRAN 
subprograms initialize the array containing T^ (EINIT) and read in the 
starting values (iNIT and RDR) . The calculation that includes thermal 
conductivity which is referred to in Chapter 3 is not discussed in detail 
here. In order to avoid undue restriction on the time step [to satisfy 
the Courant-Friedrichs-Lewy condition - see Richtmyer and Morton (196?) 
the method employed' is implicit and requires the inversion of a tridia- 
gonal matrix (dimension 8ij-6) and is rather slow. Runs with this program 
indicated thermal conduction did not change the results substantially 
so these routines are not reproduced here. 
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OE POOR quality 


TABULATION ROUTINE 


IMPLICIT REAL*8 CA-H,0~Z) 

C 

C THIS PROGRAM CALCULATES SEMI=LOGARITHMIC TABLES NEEDED FOR 

C CALCULATION OF REVERSE CURRENT HEATING WITH RESISTIVITY THAT 

C IS A FUNCTION OF THE DRIFT VELOCITY. 

C 

C THE TABLES ARE WRITTEN OUT TO LOGICAL UNIT 9 
C 

REAL*8 TL ( 820) ,TM( 820) ,VITH( 820) , CHIN (820) ,CINV( 564) 

REALMS DTO/2 431 0000000000000/, TB0/Z44 10000000000000/, 

.KAY/1 . 38054D-16/,EC/4.8029SD-10/,MH/1 . 6753089D-24/, 
.PI/Z413243F6ASSS5A31/,MP/1 . 67252D-24/, ME/9 . 109D-2S/, 

. HBAR/1 . 054 500- 27/, MU, El ON, ONE 3/ZC0555 555 55555555/, 

. C/2.997925D 10/, CSIG, LSIG,K32,CKAP,LNS1 ,LNS2, DCON, BETA, PSI 
INTEGER*4 LIM(4)/242, 242, 242, 51/ 

9001 FORMATCIOAS) 

CALCULATE CONSTANTS 

MU=(ME*MP)/(ME+MP) 

EI0N = (MU*EC^-'»''4)/(2. D0*HBAR**2) 

TI0N=(2. D0^EI0N)/(3. D0*KAY) 

CSIG = 0. 5D0'*'PI 
K32=KAY-*^DSQRT(KAY) 

CSIG=(0. 5D0TcSIG^DSQRT(CSIG^ME)*EC*EC*C)/K 32 
LSIG = (3.D0't-K32)/(2. D0^:EC*EC*EC*DSQRT(PD) 

CVITH = 2. DO^’169. DO^-KAY/MP 
GCHIN-1 . 5-+KAY 
LNS1=DL0G(LSIG) 

LNS2 = DL0G(LSIG^'4.2D5) 

DCON=1.DO/(4. 505^1. 58D5) 

CALC tables: ELECTRICAL CONDUCTIVITY (TL,TM), CRITICAL 

VELOCITY (VITH*:13.) AND INVERSE IONIZATION FRACTION (CHIN). 

DT=DT0 
T0=TB0 

DO 20 11=1,4 

K=(256*(II-1))+1 
L=K+LIM(II) 

T=T0-DT 
DO 10 J=K,L 

BETA=1 . 5SD5/T 

PSI=4.5D5/(BETA:»:DEXP(BETA)*( . 428SD0+. 5DO:<-'DLOG ( BETA) 

+ . 4698DO'<:DETA*:*^ONE3) ) 

VITH( J)=DSQRT(CVITHn) 

CHIN( J)=CCHIN:*’( 1 . DO+PSI )/PSI 
CHI=PSI/( 1 . DO+PSI ) 

IF(T.GT. 4.2D5)G0T0 5 
TM(J)=CSIG/(T:*-DSQRT(T) ) 

TL(J)=TM(J)*:(LNS1+1 .5DO:FDLOG(T)-. 5D0*DL0G(CHI ) ) 

GOTO 10 

5 TM( J)=CSIG/(T:< DSQRT(T) ) 

TL( J)=TN(J)*(LNS2+DL0G(T)-. 500*DL0G(CHI ) ) 


10 

T=T+0T 



T0 = T0:F16. 

DO 

20 

DT = DT:»:16. 

DO 


CALCULATE TABLES FOR CHI INVERSE ( ( 1+CHI )*T+CHI*TI0N AS FUNC OF T) 

TNEW=TB0-DT0 

DT=DT0 

T0=TB0 

DO 40 11=1,3 

K = (256'^(II-1)) + 1 
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L=K+LIM(II+1) 

T=TO-DT 
DO 30 J=K,L 

TCH=4. D-16*TNEW 
25 TOLD=TNEM 

BETA=1 . 58D5/T0LD 
EBETA=DEXP(BETA) 

B13=BETA**0NE3 

TEMP1=0. 4288D0+0. 5DO=*:DLOG(BETA)+. 4698D0*B13 
PSI=4. 5D5/(BETA^EBETA*TEMP1) 

C=PSI/( 1 . DO+PSI ) 

DC = DC0N*C*C'^BETA*BETA*EBETA*((1 . DO+BETA)>i; 
TEMPI + ( . 5D0-. 1566D0*B13)) 

TNEW=TOLD-( ( 1 . DO+C)-nOLD+C*TION-T)/( 1 .DO.i'C+ 
(TOLD+TION)*DC) 

IF(DA8S(TNEW-T0LD) .GT.TCH)GOTO 25 
CINV(J)-TNEW 


30 

T=T+DT 



T0=T0*16. 

DO 

40 

DT=DT*16. 

DO 


WRITE OUT TABLES 

URITE(9, 9001 )TL 

URITE(9,9001)TM . 

WRITEC9, goonviTH 

WR1TE(9,9001)CHIN 

WRITE(9, 900DCINV 

STOP 

END 
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REVERSE CURRENT HEATING MAIN ROUTINE 


IMPLICIT REAL*8 (A-H,0-Z) 

THIS PROGRAM CALCULATES REVERSE CURRENT HEATING OF A MODEL 
ATMOSPHERE READ IN AS UP TO 1024 VALUES OF TEMPERATURE (T) NUMBER 
DENSITY (N) AND DISTANCE (S) FROM THE INJECTION POINT (TOP OF 
MODEL. THE PROGRAM DOES NOT HAVE TO START AT TIME 0 (INJECTION 
TIME) AS THE CURRENT TIME, TIME STEP AND ITERATIONS TO THIS POINT 
ARE READ IN ALSO. THE PROGRAM READS IN THE MAXIMUM NUMBER OF 
ITERATIONS TO BE PERFORMED (NITER), THE ENERGETIC ELECTRON NUMBER 
FLUX (EFLUX), PSIO WHICH CORRESPONDS TO AN ENERGY CHARACTERISTIC 
OF A LOW ENERGY KNEE IN THE ENERGETIC ELECTRON DISTRIBUTION 
(SEE KNIGHT AND STURROCK 1977) , FRAC, THE MAXIMUM PERCENTAGE 
CHANGE IN THE THERMAL ENERGY CONTENT OF THE PLASMA PER HYDROGEN 
NUCLEUS ALLOWED AT ANY GRID POINT IN ONE TIME STEP, TIMMAX 
THE MAXIMUM TIME FOR THIS RUN (REAL TIME NOT COMPUTER TIME) AND 
DTOUT, THE INTERVAL AT WHICH THE ARRAYS CONTAINING THE TEMPERATURE 
AND CURRENT DENSITY AS WELL AS THE CURRENT TIME AND TIME STEP. 

CALLED SUBROUTINES: 

INIT - READS IN STARTING VALUES OF DENSITY, TEMPERATURE AND 

DISTANCE AS WELL AS TIME, TIME STEP AND NUMBER OR PREVIOUS 
ITERATIONS. 

EINIT - CALCULATES INITIAL TE DEFINED AS ( 1+CH I ) T+2*E I 0N/3*K FOR 
EACH SPATIAL GRID POINT. ENERGY INPUT INCREASES TE 
AND T IS CALCULATED FROM CHINV. 

NOUT - WRITES OUT DENSITY AND DISTANCE ARRAYS AS WELL AS INPUT 
PARAMETERS (TO FORTRAN LOGICAL UNIT 9) 

CURINT - INITIALIZATION FOR CURCAL (SEE CURCAL) 

TESTPI - INITIALIZATION FOR TESTP (SEE TESTP) 

CURCAL - CALCULATES CURRENT AS A FUNCTION OF DISTANCE USING 
STEADY STATE RESULTS OF KNIGHT AND STURROCK AND A 
RESISTIVITY THAT DEPENDS ON THE REVERSE CURRENT DRIFT 
VELOCITY . 

TESTP - UPDATES TEMPERATURE AND ADJUSTS TIME INCREMENT SO THAT 

MAXIMUM CHANGE IN TE AT ONE GRID POINT IS FRAC*TE AT THE 
GRID POINT 

DIAG - OUTPUTS A SMALL SUBSET OF THE CALCULATED CURRENT DENSITY 
AT INTERVALS DETERMINED BY VALUES IT READS FROM LOGICAL 
UNIT 5 - CAN BE RECOMPILED WITHOUT RECOMPILING THE REST 
OF THE PROGRAMS AS IT DOES NOT AFFECT CALCULATIONS. 

TOUT - WRITES OUT TEMPERATURE AND CURRENT ARRAYS AND CURRENT 
TIME, TIME STEP AND ITERATIONS TO LOGICAL UNIT 9 

CTOUT - CLOSES LOGICAL UNIT 9 (I.E. END FILE 9) 


DECLARE VARIABLES: 

REAL*8 KAY, ME, C, GAM, NORM, SFST,EFLUX, PSIO, DTOUT, 
.EL,EFACT,NEWDT 

REAL'*8 T( 1024) ,N(1024) ,J( 1024 ),S( 1024) ,OSIG( 1024) ,LN( 1024), 
.TE( 1024) , TUP ( 1024) ,SD( 1024), 

.TL(820),TM(820) ,CNV(820),VITH(S20), GINV(564) 

INTEGER:^’4 IND(2, 1024) 

5001 F0RMAT(5F7.0) 

5002 FORMAT (I 7) 
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9001 FORnATdOAS) 

c 

C INITIALIZE CONSTANTS: BOLTZMANN'S CONSTANT, ELECTRON REST MASS, 

C ELECTRON CHARGE (ESU), # ERG/KEV, SPEED OF LIGHT 

C 

KAY=1 . 380540-16 
ME=9. 10910-27 
EL=4. 802980-10 
EFACT=1 . 60210D-9 
C=2. 997925010 
C 

C READ IN MODEL PARAMETERS 
C 

READ(5,5002)NITER 

READ(5,5O01)EFLU>:,PSIO,FRAC,TIMMAX,DTOUT 
READ(8,9001)TL 
READ(8,9001)TM 
REA0(8,9001)VITH 
REA0(8,9001)CNV 
REA0(8,9001)CINV 

UNITS OF INPUT PARAMETERS ARE= 

EFLUX 1.D17 (CMT:t:2 SEC)^'»:-1, PSIO IN KEV, FRAC IN PERCENT, 

TIMMAX IN SEC (MAXIMUM TIME),DT0UT IN SEC (OUTPUT INTERVALS) 

SCALE INPUT VARIABLES 

EFLUX=EFLUX*1 .D17 
PSI0=(PSI0=«EFACT)/EL 
FRAC=FRAC=*-0.01D0 

CALCULATE CONSTANTS FOR CALCULATION OF CURRENT 

TEMP1=PSI0+PSI0 
GAM=2.5D0 
GAMM1=1,5D0 

TEMP 1=TEMP1*0SQRT( TEMPI) 

TEMP2 = PSI0'*DSQRT(PSI0) 

N0RM=(EFLUX»-GAMM1-VIETTEMP1)/EL 
TEMP1=-(EL'*-’EL‘*-N0RM)/(GAMM1=^ME*C) 

J( 1 )=TEMP1/TEMP2 
TEMP2=TEMP2*PSI0 
TEMP3=-GAM'*-TEMP1 

INITIALIZE TIME AND ITERATIONS 

NE1-!DT = 0.D0 
IITER=0 
TIM=0. DO 


INITIALIZE TEMPERATURE AND DENSITY 
NTAB=1024 

CALL INIT(T,N,S,DELT,TIM,NIT,NTAB) 

OUTIME=DTOUT+TIM 

NITER=NIT+NITER 

IITER=NIT 

T(NTAB+1)=T(NTAB) 

S(NTAB+1)=S(NTAB) + (S(NTAB)-S(NTAB-D) 

CALL EINIT(T,TE,NTAB) 

OUTPUT INPUT PARAMETERS AND INITIAL DENSITY AND DISTANCES 

CALL NOUT(EFLUX,PSIO,FRAC,TIMMAX,NTAB,N,S) 

NEED 1/N IN LOOP SO WE CHANGE N TO 1/N AND CALCULATE 
DIFFERENCES OF DISTANCES USED IN TIME STEP. 
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SFST=-S(2) 

DO 10 I=1,NTAB 
SD(I)=(S(I)-SFST)*0. 5D0 
SFST=S(I) 

LN(I)=0. 500*DL0G(N(n) 

10 N(I)=2.D0/(3.D0^KAY»^Nm) 

LN(NTAB+1 )=LN(NTAB) 

PASS ADDRESSES OF INTERPOLATION TABLES AND OTHER 
CONSTANTS TO CURCAL AND TESTP 

CALL CURINT(TEMPl,TEriP2,TEMP3,GAM,TL,TM, VITH,CNV, 

•N,LN,S0,0SIG, J,T,NTAB) 

CALL TESTPI(T,TE, J,0SIG,N, TUP, BELT, FRAC,NTAB,CINV) 

START TIME STEPPING LOOP 
1 CONTINUE 

CALCULATE CURRENT AND RESISTIVITY AT EACH GRID POINT 
CALL CURCALCOSIG, J,T,N,LN,NTAB) 

CALCULATE ONE TIME STEP WORTH OF HEATING, UPDATE TEMPERATURE 
' AND ADJUST TIME STEP ACCORDING TO FRAC 

CALL TESTP(T,TE, J,OSIG,N,TUP,DELT,FRAC,NTAB, CINV) 

WRITE OUT SOME STUFF TO MAKE SURE THINGS ARE WORKING RIGHT 

CALL DIAG(S,T,TE, J,OSIG,N,TUP,DELT,TIM,NTAB, IITER) 

STEP TIME 

'lITER=IITER+1 
IF(IITER.GT.NITER)GOTO <10 
TIM=TIM+DELT 
IFCTIM.lt. OUTIMEIGOTO 1 
0UTIME=0UTIME+DT0UT 

OUTPUT CURRENT VALUES OF TIM, TEMP AND J 

CALL TOUT (TIM,T,J, DELT, IITER, NTAB) 

IF(TIM.LT.TIMMAX)GOTO 1 
40 CALL CTOUT 

WRITECe, WHITER, DELT 

STOP 

END 

SUBROUTINE E I N I T ( T , TE , NTAB ) 

IMPLICIT REALMS (A-H,0-Z) 

REAL*8 T(NTAB) ,TE(NTAB) 

REALMS KAY/ 1 . 380 54 D- 1 6/, EC/4 . 802980- 1 0/ , MP/ 1 . 67252D-24/, 

.ME/9. 1091D-28/, 

. HDAR/1 . 05450D-27/,ONE3/ZC055555555555555/, 

. BETA, CHI 

TION=(ME=«MP)/(ME+MP) 

TION=(TION*EC**4)/(2. D0^HRAR**2) 

TI0N = (2. 00m0N)/(3. DO^-KAY) 

DO 10 1=1, NTAB 

BI:TA = 1 . 5SD5/T( I ) 

CHI=4. 5D5/(BETA^:dEXP(BETA)*( .4288D0+. 5D0^-'DL0G(BETA) 

+ . 469S'^BETA^'*0NE3) ) 

CHI=CHI/( 1 . DO+CHI) 

10 TE(I)=T(I ) + (CHI=»''(T(n+TION)) 

RETURN 

END 

SUBROUTINE I N I T ( T , N , S , D ELT , T IM, N I T , NTAB ) 

IMPLICIT REAL-t^S i:a-H,0-Z) 

REAL*S T(1024),N(1024),S(1024) 
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1001 FORMATC10A8) 

READ (10, 1001 )NTAB, NIT, DELT, TIM 
CALL RDR(T,N,S,NTAB) 

RETURN 

END 

SUBROUTINE RDR (T , N , S, NTAB ) 
REALMS T(NTAB),S(NTAB),N(NTAB) 
1001 FORMAT(IOAS) 

READ! 10, 100DN 
READCIO, 100DS 
READC 10, 100DT 
RETURN 
END 
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CURCAL 



CURCAL CSECT 
* 

* ROUGHLY EQUIVALENT TO THE FORTRAN CODE BELOW, EXCEPT 

* THE FUNCTIONS PASSED IN THE ARGUMENT LIST OF THE 

* FORTRAN ENTRY POINT CURINT (TL,TM, VTH, CNV) ARE IMPLEMENTED 

* IN LINE IN THE ASSEMBLY LANGUAGE VERSION AND A CALL TO 

* THE ASSEMBLY LANGUAGE VERSION SHOULD PASS THE ADDRESSES 

* OF SEMI-LOGARITHMIC TABLES (TL ( 820 ) , TM ( 820 ) , VTH ( 820 ) , CN V ( 820) ) 

* VIA THE ENTRY POINT CURINT RATHER THAN FUNCTION NAMES. 

* CURCAL DOES THE EQUIVALENT OF A FORTRAN RETURN 1 WHEN A 

^ VALUE OF T IS OUTSIDE THE TABULATED RANGE. THERE IS NO 

* OBVIOUS WAY TO MAKE THIS APPARENT IN THE FORTRAN VERSION. 

* THE ARGUMENTS PASSED TO CURCAL ARE IGNORED AND 

* OBTAINED FROM LOCAL STORAGE WHERE CURINT PUT THEM. 

* 

* NOTE THAT THIS MEANS CURINT MUST BE CALLED BEFORE THE 

* FIRST TIME CURCAL IS CALLED OR A REAL MESS WILL RESULT. 

* 

* 

* SUBROUTINE CURI NT (T EMP 1 , TEMP2, TEMP3, GAM, TL , TM, VTH , CNV , 

* .N,LN,SD,OSIG, J,T,NTAB) 

* 

* DECLARE VARIABLES 


* 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

* 

* 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 
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IMPLICIT REALMS (A-H,0-Z) 

REAL*8 N(NTAB) ,LN(NTAB) ,SD(NTAB) ,T(NTAB) , J(NTAB) ,OSIG(NTAB) , 
.TL,TM,VTH,CNV 

REAL*8 TEMPI, TEMP2,TEMP3, CONAN, ESU,C,CONAN1,B,CHT, 

. VC,VITH,TSI , JA,NCON,NCON1,ME,MI,PI, JCON, JC1 
DATA C/2. 997925D10/, ESU/4 . 80298D- 1 0/, ME/9 . 109D-28/, 

.MI/1 . 67252D-24/,PI/Z413243F6A8885A30/, ERR/1 . D-6/ 

INITIALIZE CONSTANTS FOR CURCAL 

INTEGER*4 COUNT/30/ 

B=1.D2 

C0NAN = 6. D-2'«DSQRT(ME/MI )*B/(4.D0*PI*ESU) 

TCON=( (GAM- 1 . DO)/GAM) 

C0NAN1=-TC0N^TEMP3^C0NAN^:ESU/C 

TC0N=TC0N'FTEMP1 

CE=-C/ESU 

RETURN 

ENTRY CURCALCOSIG, J,T,N,LN,NTAB,*) 

CALCULATE CURRENT AND RESISTIVITY (OSIG) FOR FIRST POINT 

0SIG(1)=TL(T( 1 ) )+TM(T( 1))*LN(1) 

CHT=CNV(T( 1)) 

VITH=VTH(T( 1 ) ) 

VD=J( 1 )^'N( 1 l^'CHT-^CE 
IFCVD.lt. VITHIGOTO 10 

THE NEXT STATEMENT CALCULATES THE ANOMALOUS PART OF THE 
RESISTIVITY IF THE DRIFT VELOCITY IS GREATER THAN THE 
CRITICAL VELOCITY FOR THE ONSET OF TURBULENCE 

0SIG(1)=0SIG( 1)+C0NAN*N(1)*CHT*(1 .DO-VITH/VD) 

INITIALIZE TSI TO ZERO 

10 TSI=O.DO 

DO 20 I=2,NTAB 

OSIG(I)=TL(T(I))+TM(T(I))*LN(I) 

TSI=TSI + (0SIG(I-1)+0SIG(I))*SD(I ) 

JC0N = (TEMP2+TEMP3*TSI ) 
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* F35C1=1 .DO/JCON 

* F35C2=F35(F35C1) 

* CHT = CNV(T(D) 

* VITH=VTH(T(D) 

* J(I)=TEI1P1*F35C2 

* NCONO=NCI)*CHT 

* CJA=NCONO-^SD(I)*CONAN 

* NCON=NCONO't-’CE 

* VD=JT>^NCON 

* IF(VD.LT.VITH)oOTO 20 

* THIS SECTION CALCULATES ANOMALOUS PART OF RESISTIVITY USING 

* NEUTON'S METHOD FOR INTERIOR GRID POINTS - SKIP IF 

* DRIFT VELOCITY IS LESS THAN CRITICAL VELOCITY 

* 

* JA=(12. DO/13. 00)*VITH/NC0N 

^ JC0N=JC0N+CJA*TEMP3 

* CJA=CJA*JA 
NCON1=NCON*CONAN1 

* J(I)=JA'«(1.D0+(J(n~JA)/(NC0N1*SD(I)*F35C1*J(I)+JA)) 

* DO 15 K=1 -COUNT 

* CJAT = CJA/J(n 

* TC1-TCON-CJAT 

* JC1=JCON-CJAT 

* JT=JCI)^'( (TC1+TEMP1*'JCn/(J(I)*JCl*F35(JC1)+TC1)) 

* IF(DABS((JT-J(in/JT).LE.ERR)60T0 16 

* 15 J(I)=JT 

^ 16 TSIG=NCONO'<CONAN^(1.DO-JA/J(I)) 

* TSI=TSI+TSIG^SD(I) 

* OSIG(n=OSIG(I)+TSIG 

* 20 CONTINUE 

* RETURN 

^ END 

* 



USING 

=^,15 



B 

CFIRST 

BRANCH AROUND NAME, OTHER ENTRY ETC 


DC 

X'06' 



DC 

CL7'CURCAL ' 



ENTRY 

CURINT 



USING 

*,15 


CURINT 

B 

IFIRST 

BRANCH AROUND NAME, SAVE AREA ETC. 


DC 

X'06' 



DC 

CL7'CURINT ' 


AREA 

DS 

18F 


REG1 

DC 

AL-KARGA) 

R1 (ADDR ARG LIST) 

REG2 

DC 

F ' 0 ' 

R2 BASE T 

REG3 

DC 

F'O' 

R3 (BASE OSIG - 8 - BASE T) 

REG4 

DC 

F 'O' 

RA (BASE J - BASE T) 

REG5 

DC 

F'O' 

R5 (BASE N - BASE T) 

REG6 

DC 

F'O' 

R6 (BASE LN - BASE T) 

REG7 

DC 

F'O' 

R7 (BASE SO - BASE T) 

EIGHT 

DC 

F'8' 

R8 (INCREMENT - 8) 

REG9 

DC 

F'O' 

R9 (BASE T + 8*(NTAB-1) COMPARAND) 

TLA 

DC 

F'O' 

R10 (BASE TL TABLE CHANGES IN LOOP) 

TMA 

DC 

F'O' 

R11 (BASE TM TABLE CHANGES IN LOOP) 

VTHA 

DC 

F'O' 

BASE VTH TABLE 

CNVA 

DC . 

F'O' 

BASE CNV TABLE 

COUNT 

DC 

F '30' 

MAX # OF NEWTON'S METHOD ITERATIONS 

ARGA 

DC 

X'SO' 

ARGUMENT LIST (ONE LONG) 


DC 

AL3CF35A) 

ADDR ARGUMENT 

TBND 

DC 

X'00000330' 


TDISP 

DC 

X'0000‘1410' 



CHOP 

0, 8 

FORCE DOUBLE WORD ALIGNMENT 

WM1 

DC 

D ' 0. ' 


WP1 

DC 

D '0. ' 


NM3 

DC 

D ' 0. ' 


WPS 

DC 

D'O. ' 


FLOAT 

DC 

X'4000000000000000' 

CONAN 

DC 

D'O.' 
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CONANl 

TEMPI 

TEMP2 

TEMP3 

GAM 

NCON 

NCONO 

ICON 

JCON 

CJA 

TCI 

JA 

TSI 

ERR 

C 

ESU 

ME 

M! 

B 

PI 

CE 

F35C1 
F35C2 
F35A 
I FIRST 





DC 

poo;? W 

D'O. ■ ^ 

DC 

D'O.' 


DC 

D'O. ' 


DC 

D'O.' 


DC 

D'O. ' 


DC 

D'O.' 


DC 

D'O.' 


DC 

D'O.' 


DC 

D'O.' 


DC 

D 'G. ' 


DC 

D'O.' 


DC 

D'O.' 


DC 

D'O.' 


DC 

0 ' 1 .E-6' 

ERROR TOLERANCE FOR NEWTON 

DC 

D '2.997925E10' 


DC 

D'4.80298E-10' 


DC 

D '9. 1091E-28’ 


DC 

D'1.67252E-24' 


DC 

O' 1.E2' 


DC 

X'‘4132‘13F6A8885A30' 

DC 

D'O.' 


DC 

D'O. ' 


DC 

D'O. ' 


DC 

D'O.' 


STM 

14, 12,12(13) 

SAVE CALLING ROUTINE'S GPR 

LR 

2,13 

R2 <= ADDR OLD SAVE AREA 

LA 

13, AREA 

R13 <= ADDR NEW SAVE AREA 

DROP 

15 

R15 NO LONGER BASE REG 

USING 

AREA, 13 

R13 NEW BASE REG 

ST 

2.4(13) 

LINK SAVE AREAS 

ST 

13,8(2) 


LM 

2, 9, 0(1) 

R2-R9 <= ADDR'S 1ST 8 ARG ' 

LD 

0,0(5) 

FO <= GAM 

ST 

6, TLA 

TLA <= BASE TL TABLE 

LDR 

2,0 

F2 <= GAM 

ST 

7,TMA 

TMA <= BASE TM TABLE 

LD 

4,0(4) 

F4 <= TEMP 3 

ST 

8, VTHA 

VTHA <= BASE VTH TABLE 

STD 

4.TEMP3 

TEMP3 (LOCAL) <= TEMP3 

ST 

9, CNVA 

CNVA <- BASE CNV TABLE 

LD 

6,0 (3) 

F6 <= TEMP2 

LM 

4,10,32(1) R4-RT0 ADDR'S REST OF ARC'S 

LD 

4,0(2) 

F4 <= TEMPI 

L 

10,0(10) 

RIO <= NTAB 

STD 

6, TEMP2 

TEMP2 (LOCAL) <= TEMP2 

SLA 

10,3 

RIO <= NTAB*8 

STD 

4, TEMPI 

TEMPI (LOCAL) <= TEMPI 

S 

10. EIGHT 

RIO <= 8'K(NTAB-1) 

SD 

2,=D' 1. ' 

F 2 < = GAM- 1 . D 0 

SR 

4.9 

R4 <= BASE N - BASE T 

DDR 

2,0 

F2 <= (GAM-1 . DO)/GAM 

SR 

5,9 

R5 <= BASE LN - BASE T 

STD 

2, GAM 

GAM <= (GAM-1 .00)/GAM 

SR 

6,9 

R6 <= BASE SD - BASE T 

LD 

0,ME 

FO <= ME 

L 

1,REG1 

R1 <= ADDR ARG LIST 

DD 

0,MI 

FO <= ME/MI 

L 

15, =V(DSQRT) 

R15 <= ENTRY ADDR DSQRT 

STD 

0, F35A 

F35A <= ME/MI 

BALR 

14,15 

FO <= DSQRT (ME/MI) 

SR 

8,9 

R8 <= BASE J - BASE T 

MD 

0, B 

FO <= DSQRT (ME/MI ) •'*'6 

SR 

7, 9 

R7 <= BASE OSIG - BASE T 

MD 

0,=D'6.E-2' 

FO <= 6.D-2*DSQRT(ME/MI)*B 

AR 

10,9 

RIO <= BASE T + S:^’(NTAB-1) 

LD 

2,=D'4. ' 

F2 <- 4. DO 

MO 

2, PI 

F2 <= 4.D0'*-'PI 

S 

7, EIGHT 

R7 <= BASE OSIG - 8 - BASE 

MD 

2, ESU 

F2 <= 4.D0*PI*ESU 


'S METHOD 
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ST 

9,REG2 

REG2 <= BASE T 


DDR 

0,2 

F2 <= 6.D-2*DSQRT(ME/MI)*B/(4.D0*PI*ESU) 


ST 

7.REG3 

REG3 <= BASE OSIG - 8 - BASE T 


STD 

0, CONAN 

CONAN <= 6.D-2*DSQRT(ME/MI)*B/(4.00*PI*ESU) 


ST 

8,REG4 

REG4 <= BASE J - BASE T 


LD 

4, GAM 

F4 <= ( (GAM-1 .DO)/GAM) 


MD 

0,TEMP3 

FO <= TEMP3*C0NAN 


ST 

4, REGS 

REGS <= BASE N - BASE T 


MDR 

0,4 

FO <= ((GAM-1. D0)/GAM)*TEMP3*C0NAN 


ST 

5,REG6 

REG6 <= BASE LN - BASE T 


LD 

2,C 

F2 <= C 


ST 

6,REG7 

REG7 <= BASE SD - BASE T 


LCDR 

2,2 

F2 <= -C 


DDR 

0,2 

FO <= -((GAM-1. D0)/GAM)*TEMP3*C0NAN/C 


ST 

10,REG9 

REG9 <= BASE T + 8*(NTAB-1) COMPARAND 


MD 

4, TEMPI 

F4 <= TEMP1*((GAM-1)/GAM) = TCON 


DO 

2,ESU 

F2 <= -C/ESU 


STD 

4,TC0N 

TCON <= TEMP1*((GAM-1)/GAM) 


L 

10,4(13) 

RIO <= ADDR OLD SAVE AREA 


MD 

0,ESU 

FO <= -((GAM-1. D0)/GAM)*TEMP3*C0NAN*ESU/C 


STD 

2,CE 

CE <= -C/ESU 


LM 

14,10,12(10) 

GPR'S RESTORED 


STD 

0, C0NAN1 

C0NAN1 <= ((GAM-1. D0)/GAM)*TEMP3*C0NAN 


L 

13,4(13) 

R13 <= ADDR OLD SAVE AREA 


SR 

15,15 

R15 <= 0 (RETURN CODE) 


MV I 

12( 13) ,X 'FF ' 

INDICATE CONTROL RETURNED 


BR 

14 

RETURN 


DROP 

13 



USING 

CURCAL, 15 


CFIRST 

STM 

14,12,12(13) 

SAVE CALLING ROUTINE'S GPR'S 


LR 

2, 13 

R2 <= ADDR OLD SAVE AREA 


LA 

13, AREA 

R13 <= ADDR NEW SAVE AREA 


DROP 

15 

R15 NO LONGER BASE REG 


USING 

AREA, 13 

R13 NEW BASE REG 


ST 

2,4(13) 

LINK SAVE AREAS 


ST 

13,8(2) 



LM 

1 , 1 1 ,REG1 

SET UP GPR'S 


LH 

12,0(2) 

R12 <= HIGH ORDER 2 BYTES OF T(1) 


MVC 

FLOAT+1 (6) , 2(2) FLOAT <= FRACTIONAL DISPLACEMENT 


S 

12,TDISP 

REDUCE R12 BY TDISP. NOW # DOUBLE WORDS 




FROM BASE OF INTERPOLATION TABLES 


LD 

4, FLOAT 

F4 <= FRACTION 0 LE FRAC LT 1 


BM 

BADT 

IF RESULT NEGATIVE - OUT OF RANGE 




GOTO BADT 


C 

12, TEND 

IF R12 GREATER THAN TBND 


BH 

BADT 

GOTO BADT 


SO 

4,=D' .5' 



SLA 

12,3 

R12 <= R2*8 NOW BYTE DISPLACEMENT 




FROM BASE OF INTERPOLATION TABLES 


NOW COMPUTE WEIGHTS 

FOR CUBIC INTERPOALTION OF 


FUNCTIONS OF T 



LDR 

2,4 

F2 <= X 


MDR 

4*4 

F4 <= X**2 = X2 


HDR 

4,4 

F4 <= X2/2 


SD 

4,=D' 1. 125' 

F4 <= X2/2 - 9/8 


LDR 

6,4 

F6 <= X2/2 - 9/8 


HDR 

4,4 

F4 <= X2/4 - 9/16 


MDR 

6,2 

F6 <= X3/2 - 9X/8 


LCDR 

0,4 

FO <= -X2/4 + 9/16 


ADR 

0,6 

FO <= X3/2 - X2/4 - 9X/8 + 9/16 


STD 

0,WM1 

WM1 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 




ING TO CLOSEST SMALLER VALUE OF T. 


LCDR 

0,4 

FO <= -X2/4 + 9/16 


SDR 

0,6 

FO <= -X3/2 - X2/4 + 9X/8 + 9/16 


STD 

0,WP1 

WP1 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 




ING TO CLOSEST LARGER VALUE OF T. 


ADR 

6,2 

F6 <= X3/2 - X/8 
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OB exjn Quyjjijf 


AD 

4,=D ’ . 5' 

F4 <= X2/4 - 1/16 


MD 

6, =X '4055555555555555' F6 <= X3/6 - X/24 

LDR 

0,4 

FO <= X2/4 - 1/16 


SDR 

0,6 

FO <= -X3/6 + X2/4 

+ X/24 - 1/16 

ADR 

6 , 4 

F6 <= X3/6 + X2/4 

- X/24 - 1/16 

STD 

0,14M3 

UM3 <= HEIGHT FOR 

SMALLEST VALUE UF T 

STD 

6,WP3 

UP3 <= HEIGHT FOR 

LARGEST VALUE OF T 


* NOW CALCULATE INTERPOLATED VALUES OF TL AND TM 

* (HAVE UEIGHTS FOR TABLES ENTRIES 1 C 4 IN FPR'S 0 £ 6) 


LOR 

4,0 

F4 

<= HEIGHT 1 

LDR 

2,6 

F2 

<= HEIGHT -4 

MD 

0,0(10, 12) 

FO 

<= HEIGHT 1 * TL 1 

MD 

2,24(10, 12) 

F2 

<= HEIGHT 4 * TL 4 

MD 

4,0(11, 12) 

F4 

<= HEIGHT 1 * TM 1 

MD 

6,24(11,12) 

F6 

<= HEIGHT 4 * TM 4 

ADR 

0,2 

FO 

<= HI^TLI + H4»:tL4 

ADR 

4,6 

F4 

<= H1"TM1 + H4'«TM4 

LD 

2,HM1 

F2 

<= HEIGHT 2 

LDR 

6,2 

F6 

<= HEIGHT 2 

MD 

2,8(10, 12) 

F2 

<= H2nL2 

MD 

6,8(11, 12) 

F6 

<= H2*TM2 

ADR 

0,2 

FO 

<= HinU + H4*TL4 + H2^TL2 

ADR 

4,6 

F4 

<= HinMl + H4*TM4 + H2*TM2 

LD 

2,HP1 

F2 

<= HEIGHT 3 

LDR 

6,2 

F6 

<= HEIGHT 3 

MD 

2, 16(10, 12) 

F2 

<= H3nL3 

MD 

6,16(11,12) 

F6 

<= H3nM3 

ADR 

0,2 

FO 

<C= INTERPOLATED VALUE OF TL 

ADR 

4,6 

F4 

<= INTERPOLATED VALUE OF TM 

MD 

4, 0(6, 2) 

F4 

<= TM*LN(1) 

ADR 

4,0 

F4 

<= TL + TM^Nd) 

L 

10, CNVA 

RIO 

1 <= BASE ADDR CHINV TABLE 

STD 

4, 8(3, 2) 

OSIG(l) <= TL + TMTLN(I) 

L 

1 1 , VTHA 

R11 

<= BASE VTH TABLE 

LD 

0,HM3 

FO 

<= HEIGHT FOR SMALLEST VALUE OF T 

LD 

6,WP3 

F6 

<= WEIGHT FOR LARGEST VALUE OF T 

NOH CALCULATE INTERPOLATED VALUES OF CH AND VT 

(HAVE HEIGHTS 

; FOR TABLES 

ENTRIES 1 £ 4 IN FPR'S 0 £ 6) 

LDR 

4,0 

F4 

<= HEIGHT 1 

LDR 

2,6 

F2 

<= WEIGHT 4 

MD 

0,0(10, 12) 

FO 

<= HEIGHT 1 CH 1 

MD 

2,24(10,12) 

F2 

<= HEIGHT 4 * CH 4 

MD 

4,0(11. 12) 

F4 

<= HEIGHT 1 * VT 1 

MD 

6, 24(11, 12) 

F6 

<= WEIGHT 4 * VT 4 

ADR 

0,2 

FO 

<= WI^'CHI + W4*CH4 

ADR 

4, 6 

F4 

<= W1TVT1 + U4*VT4 

LD 

2,UM1 

F2 

<= HEIGHT 2 

LDR 

6,2 

F6 

<= HEIGHT 2 

MO 

2,8(10,12) 

F2 

<= H2'»-CH2 

MD 

6,8(11,12) 

F6 

<= U2'*^VT2 

ADR 

0,2 

FO 

<= U1*CH1 + U4'^-CH4 + H2*CH2 

ADR 

4,6 

F4 

<= H1*VT1 + W4'»^VT4 + W2'^VT2 

LD 

2,WP1 

F2 

<= WEIGHT 3 

LOR 

6,2 

F6 

<= WEIGHT 3 

MD 

2, 16(10,12) 

F2 

<= H3TCH3 

MD 

6, 16(11, 12) 

F6 

<= U3t-VT3 

ADR 

0,2 

FO 

<= INTERPOLATED VALUE OF CH 

ADR 

4, 6 

F4 

<= INTERPOLATED VALUE OF VT 

LDR 

2,0 

F2 

<= CHT 

SDR 

6 , 6 

F6 

<= 0. DO 

MD 

2, CE 

F2 

<= chtt-ce 

L 

10, TLA 

RIO 

! <= BASE TL TABLE 

MD 

2, 0(5,2) 

F2 

<= N(1)'*-CHT^CE 

L 

1 1,TMA 

RH 

<= BASE TM TABLE 

STD 

6,TSI 

TSI 

<= O.DO 
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MD 

2, 0(4, 2) 

F2 <= J(1)*N(1)*CHT*CE = VD 


COR 

2.4 

IF(VO. LT. VITH) 


BL 

ARND 

GOTO ARND 


LD 

6, CONAN 

F6 <= CONAN 


MDR 

6,0 

F6 <= CONAN^CHT 


MD 

6, 0(5, 2) 

F6 <= C0NAN*N(1)*CHT 


LD 

0,=D‘ 1. ' 

FO «C= 1 . DO 


DDR 

4,2 

F4 <= VITH/VD 


SDR 

0,4 

FO <= 1 . DO-VITH/VD 


MDR 

6,0 

F6 <= C0NAN*N(1)^CHT»(1. DO-VITH/VD) 


AD 

6, 8(3, 2) 

F6 <= OSIG( 1)+C0NAN*N(1)*CHT*( 1 .DO-VITH/VD) 


STD 

6, 8(3,2) 

0SIG(1X=0SIG( 1 )+CONAN'«N( 1)*CHT^'( 1 . DO-VITH/VD: 

ARND 

AR 

2,8 

R2 <= ADDR T(2) 

LOOP 

LH 

12,0(2) 

R12 <= HIGH ORDER 2 BYTES OF T(D 


MVC 

FL0AT+U6) 

,2(2) FLOAT <= FRACTIONAL DISPLACEMENT 


S 

12,TDISP 

REDUCE R12 BY TDISP. NOW # DOUBLE WORDS 




FROM BASE OF INTERPOLATION TABLES 


LD 

4, FLOAT 

F4 <= FRACTION 0 LE FRAC LT 1 


BM 

BADT 

IF RESULT NEGATIVE - OUT OF RANGE 




GOTO BADT 


C 

12,TBND 

IF R12 GREATER THAN TBND 


BH 

BADT 

GOTO BADT 


SD 

4,=D ' .5' 

F4 <= FRAC -.5 LE FRAC LT .5 


, SLA 

12, 3 

R12 <= R2*8 NOW BYTE DISPLACEMENT 




FROM BASE OF INTERPOLATION TABLES 


NOW 

COMPUTE WEIGHTS FOR CUBIC I NTERPOALT I ON OF 


FUNCTIONS OF T 



LDR 

2,4 

F2 <= X 


MDR 

4,4 

F4 <= X=*:*2 = X2 


HDR 

4,4 

F4 <= X2/2 


SD 

4,=0 ' 1 . 125 

' F4 <= X2/2 - 9/8 


LDR 

6,4 

F6 <= X2/2 - 9/8 


HDR 

4,4 

F4 <= X2/4 - 9/16 


MDR 

6,2 

F6 <= X3/2 - 9X/8 


LCDR 

0,4 

FO <= -X2/4 + 9/16 


ADR 

0,6 

FO <= X3/2 - X2/4 - 9X/8 + 9/16 


STD 

0.WM1 

WM1 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 

:4c 



ING TO CLOSEST SMALLER VALUE OF T. 


LCDR 

0,4 ■ 

FO <= -X2/4 + 9/16 


SDR 

0,6 

FO <= -X3/2 - X2/4 ->■ 9X/8 + 9/16 


STD 

0.WP1 

WP1 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 




ING TO CLOSEST LARGER VALUE OF T. 


ADR 

6,2 

F6 <= X3/2 - X/8 


AO 

4,=D' .5' 

F4 <= X2/4 - 1/16 


MD 

6,=X’4055555555555555' F6 <= X3/6 - X/24 


LDR 

0,4 

FO <= X2/4 - 1/16 


SDR 

0, 6 

FO <= -X3/6 + X2/4 + X/24 - 1/16 


ADR 

6,4 

F6 <= X3/6 + X2/4 - X/24 - 1/16 


STD 

0,WM3 

WM3 <= WEIGHT FOR SMALLEST VALUE OF T 


STD 

6,WP3 

WP3 <= WEIGHT FOR LARGEST VALUE OF T 


NOW CALCULATE INTERPOLATED VALUES OF TL AND TM 


WEIGHTS 

FOR TABLES 

ENTRIES 

18 4 

IN 

FPR'S 08 6) 

LDR 

4,0 

F4 

< = 

WEIGHT 

1 


LDR 

2,6 

F2 

< = 

WEIGHT 

4 


MD 

0,0(10, 12) 

FO 

<- 

WEIGHT 

1 

* TL 1 

MD 

2,24(10,12) 

F2 

< = 

WEIGHT 

4 

TL 4 

MD 

4,0(11,12) 

F4 

< = 

WEIGHT 

1 

* TM 1 

MD 

6,24(11, 12) 

F6 

< = 

WEIGHT 

4 

* TM 4 

ADR 

0,2 

FO 

< = 

Wl't-TLI 

+ 

W4*TL4 

ADR 

4,6 

F4 

< = 

WI-^TMI 

+ 

W4*TM4 

LD 

2,WM1 

F2 

< = 

WEIGHT 

2 


LDR 

6,2 

F6 

< = 

WEIGHT 

2 


MD 

2,8(10,12) 

F2 

< = 

W2TTL2 



MD 

6,8(11,12) 

F6 

< = 

W2*TM2 



ADR 

0,2 

FO 

< = 

wmLi 

+ 

W4*TL4 + W2*TL2 
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OB EooR dUAury 


ADR 

4.6 

F4 <= W1*TM1 + W4*TM4 + W2*TM2 

LD 

2.WP1 

F2 <= WEIGHT 3 

LDR 

6,2 

F6 <= WEIGHT 3 

MD 

2,16(10,12) 

F2 <= W3=PTL3 

«D 

6, 16(11,12) 

F6 <= W3nM3 

ADR 

0,2 

FO <= INTERPOLATED VALUE OF TL 

ADR 

4,6 

F4 <= INTERPOLATED VALUE OF TM 

MO 

4, 0(6, 2) 

F4 <= TM*LNCI) 

ADR 

0,4 

FO <= TL + TM'-i-LNm 

L 

10, CNVA 

RIO <= BASE ADDR CHINV TABLE 

STO 

0,8(3, 2) 

OSIG(I) <= TL + TM-t^LNd) 

L 

1 1 , VTHA 

R1 1 <= BASE VTH TABLE 

AO 

0,0(3, 2) 

FO <= OSIG(I)+OSIG(I-1) 

MD 

0,0(7, 2) 

FO <= (OSIG(I)+OSIG(I-1))^'SO(I) 

AD 

0,TSI 

FO <= TSI (0SIG(I)+nsiG(I-1))*SD(I) 

STO 

0,TSI 

TSI <= TSI + ai5IG(i)+0SIG(I-l))^S0(I) 

MD 

0,TEMP3 

FO <= TSI*TEMP3 

AD 

0,TEMP2 

FO <= TEMP2 + TEMP3*TSI 

STD 

0, JCON 

JCON <= TEMP2+TEMP3^’TSI 

LD 

6,=D' 1 . ' 

F6 <= 1.00 

DDR 

6,0 

F6 <= 1.D0/(TEMP2+TEMP3*TSI) 

L 

15,=V(F35) 

R15 <= ENTRY ADDR F35 

STO 

6,F35C1 

F35C1 <= 1.00/(TEMP2+TEf’lP3'«TSI) 

STD 

6, F35A 

F35 ARGUMENT <= 1 . D0/(TEMP2+TEMP3nsI ) 

BALR 

14,15 

FO <= F35(T.D0/CTEMP2+TEMP3nsn) 

STD 

0, F35C2 

F35C2 <= F35C1 .D0/(TEMP2+TEMP3-TSI)) 

LD 

0,WM3 

FO <= WEIGHT FOR SMALLEST VALUE OF T 

LD 

6,UP3 

F6 <= WEIGHT FOR LARGEST VALUE OF T 


* NOW CALCULATE INTERPOLATED VALUES OF CH AND VT 

* (HAVE WEIGHTS FOR TABLES ENTRIES IE 4 IN FPR'S 0 £ 6) 


* 

LDR 

4,0 

F4 <= WEIGHT 1 


LDR 

2,6 

F2 <= WEIGHT 4 


MD 

0, 0( 10, 12) 

FO <= WEIGHT 1 * CH 1 


MD 

2,24(10.12) 

F2 <= WEIGHT 4 * CH 4 


MD 

4,0(11.12) 

F4 <= WEIGHT 1 * VT 1 


MO 

6, 24(11,12) 

F6 <= WEIGHT 4 •'*: VT 4 


ADR 

0,2 

FO <= WI-'-'CHl + W4*CH4 


ADR 

4,6 

F4 <= Wl't-VTI + W4»^VT4 


LD 

2.WM1 

F2 <= WEIGHT 2 


LDR 

6,2 

F6 <= WEIGHT 2 


MD 

2,8(10, 12) 

F2 <= W2»-CH2 


MD 

6,8(11,12) 

F6 <= W2^VT2 


ADR 

0,2 

FO <= Wr+CH1 + W4t-'CH4 + W2*CH2 

ADR 

4,6 

F4 <= WdVTI + W4S-VT4 + W2*VT2 

LD 

2.WP1 

F2 <= WEIGHT 3 


LDR 

6.2 

F6 <= WEIGHT 3 


MD 

2, 16(10,12) 

F2 <= W3'CH3 


MD 

6,16(11,12) 

F6 <= W3^VT3 


ADR 

0,2 

FO <= INTERPOLATED VALUE OF 

CH 

ADR 

4,6 

F4 <= INTERPOLATED VALUE OF 

VT 

MD 

0, 0(5,2) 

FO <= N(!)*Ch'T 


LDR 

6,0 

F6 <= N(I)*CHT = NCONO 


STO 

0, NCONO 

NCONO <= N(I)'FCHr 


MD 

6, 0(7, 2) 

F6 <= NCONO-t-SOd) 


MD 

0, CE 

FO <= Nd):t-'CHT*CE 


MO 

6, CONAN 

F6 <= CONAN 'NCONO-'^SD(I) 


STD 

0, NCON 

NCON <= Nd)*CHT^CE 


STD 

6, CJA 

CJA <= CONAN't-HCONO^-SDCI) 


LD 

2, TEMPI 

F2 <= TEMPI 


MD 

2, F35C2 

F2 <= TFMPr<T35C2 = J 


L 

1 1 , TMA 

R n <= BASE TM TABLE 


STD 

2, 0(4, 2) 

Jd) <= TEMP1^T35C2 


MDR 

2,0 

F2 <= JT'-'-NCON = VD 


L 

10, TLA 

RIO <= BASE TL TABLE 


CDR 

2,4 

IF(VD.LT.VITH) 


BL 

* 

ARND1 

GOTO ARNOl 
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* 

* 

* 

* 

* 


NEWT 


* 




OUT 


* 


THE FOLLOWING CALCULATES THE CURRENT AND RESISTIVITY IN THE 
CASE THE DRIFT VELOCITY EXCEEDS 13 TINES THE ION THERMAL 

VELOCITY IN THIS CASE AT LEAST PART OF THE SLAB 

IS CHARACTERIZED BY ANOMALOUS RESISTIVITY 


LD 

2, CJA 

F2 <= CJA 

MO 

2,TEMP3 

F2 <= TEMP3*CJA 

LDR 

6,4 

F6 <= VITH 

STD 

2, CJA 

CJA <= TEMP3*CJA 

AD 

2, JCON 

F2 <= JCON + TENP3*CJA 

STD 

2, JCON 

JCON <= JCON + TEMP3'*CJA 

DDR 

6,0 

F6 <= ( 12./13. )*VITH/NCON = JA 

LD 

2, CJA 

F2 <= CJA 

MO 

0, C0NAN1 

FO <= NC0N*C0NAN1 = NC0N1 

MDR 

2,6 

F2 <= CJA*JA 

STD 

6 , JA 

JA <= (12./13. )*VITH/NCON 

STD 

2, CJA 

CJA <= CJA*JA 

LD 

2,0(4,21 

F2 <= J(I) 

MD 

0,0(7,21 

FO <= NC0N1*S0(I) 

MD 

0,F35C1 

FO <= NC0N1*SD(I)*F35C1 

MDR 

0,2 

FO <= NC0N1^SD(I)*F35C1*J(I) 

ADR 

0,6 

FO <= NC0N-'^SD(I)*F35C1*J(I)+JA 

SDR 

2,6 

F2 <= J(I)-JA 

DDR 

2,0 

F2 <= (J(n-JA)/(NC0N1*SD(I)*F35C1*J(I)+JA) 

AD 

2, =D ' 1 . ' . 

F2 <= 1.D0 + (J(I)-JA) 

/(NC0Nl*SD(I)*F35C1*J(n+JA) = MULCON 

MDR 

2,6 

F2 <= JA^MULCON 

STD 

2, 0(4, 2) 

J(I) <= JA'^'MULCON 

L 

1 1 , COUNT 

R11 <= COUNT(MAX # NEWTON'S METHOD STEPS) 

LD 

6, CJA 

F6 <= CJA 

DDR 

6,2 

F6 <= CJA/Jd) = CJAT 

L 

1S,=V(F35) 

R15 <= ENTRY ADDRESS F35 

LD 

0, JCON 

TO <= JCON 

SDR 

0,6 

FO <= JCON - CJAT = JC1 

STD 

0, F35A 

F35A <= JC1 

no 

6,TC0N 

F6 <= TCON'*:CJAT = TCI 

STD 

6,TC1 

TCI <- TCON^CJAT 

BALR 

14,15 

FO <= F35(JC1) 

LD 

4, TCI 

F4 <= TCI 

LD 

6,F35A 

F6 <= JC1 

MDR 

0,6 • 

FO <= F35(JC1)*JC1 

LD 

2, 0(4, 2) 

F2 <= J(I) 

MDR 

0,2 

FO <= J(I)*JC1*F35(JC1) 

MD 

6, TEMPI 

F6 <= TEMPI^JCI 

ADR 

6,4 

F6 <= TCI + TEMPI^JCI 

ADR 

0,4 

FO <= TC1+J(I)^JC1*F35(JC1) 

DDR 

6,0 

F6 <= (TC1+TEMP1*JC1)/(TC1+J(I)*JC1*F35(JCD: 

MDR 

2,6 

F2 <= J(I)*(TC1 + TEMP1't-JC1)/ 
(TCl+J(I)=tJCl*F35(JC1)) 

LDR 

4,2 

F4 <= NEW ESTIMATE OF J(I) 

LPOR 

0,2 

FO <= DABS(NEW ESTIMATE OF J(D) 

SD 

4, 0(4, 2) 

F4 <= NEW ESTIMATE - OLD ESTIMATE 

MO 

0, ERR 

FO <= ERR^'DABS(NEW ESTIMATE) 

LPOR 

4,4 

F4 <= DABS(NEW ESTIMATE - OLD ESTIMATE) 

CDR 

0,4 

IF(OABS(NEW ESTIMATE - OLD ESTIMATE) 
/DABS(NEW ESTIMATE) . LE . ERR) 

STD 

2,0(4,21 

J(I) <= J(I)*(TC1+TEMP1*JC1)/ 
(TC1+J(I)'*:JC1*F35(JC1)) 

BNL 

OUT 

GOTO OUT 

BCT 

1 1,NEWT 


LD 

4, JA 

F4 <= JA 

DDR 

4,2 

F4 <= JA/J(I) 

L 

1 1 ,TMA 

R1 1 <= BASE ADDR TM ARRAY 

LD 

2,=D ' 1 . ' 

F2 <= 1.D0 

SDR 

2,4 

F2 <= 1.D0 - VITH/VD 

MD 

2,NC0N0 

F2 <= NCONO*(1,DO-VITH/(J(I)=«NCON)) 

MD 

2, CONAN 

F2 <= CONAN*NCONO*(1.DO-VITH/(J(I)*NCON)) 
= TSIG 

LDR 

4,2 

F4 <= TSIG 
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ARN01 


BADT 


AD 

2, 8(3, 2) 


F2 <= OSIG(I) + TSIG 

MD 

4, 0(7,2) 


F4 <= TSIG^SDd) 

AD 

4,TSI 


F4 <= TSI + TSIG=*'SD(I) 

STD 

2, 8(3, 2) 


OSIG(I) <= OSIGd ) V TSIG 

STD 

4,TSI 


TSI <= TSI + TSIG^SOd) 

BXLE 

2,8, LOOP 


I <= 1+1 AND GOTO LOuP IF NOT DONE 

L 

13,4(13) 


R13 <= ADDR OLD SAVE AREA 

LM 

14,12, 12(13) 

GPR'S RESTORED 

SR 

15,15 


R15 <= 0 (RETURN CODE) 

MVI 

12(13), X 

'FF' 

INDICATE CONTROL RETURNED 

BR 

14 


RETURN 

L 

13,4(13) 


R13 <= ADDR OLD SAVE AREA 

LM 

14, 12,12(13) 

GPR'S RESTORED 

LA 

15,4 


R15 <= 4 (RETURN CODE) 

MVI 

12( 13),X ‘ 

'FF' 

INDICATE CONTROL RETURNED 

BR 

14 


RETURN 

END 





DKiGINXn ESGE I3f 
RE POOR QUAUTX 
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F35 


F35 

* 

* 

* 

* 

* 

* 

. 1 : 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

*. 


AREA 

FIRST 


NONNEG 


* 


NOEXTD 


CSECT 

REAL FUNCTION F35*8(X) 

F35 RETURNS THE 3./S. POWER OF THE ARGUMENT IF 

THE ARGUMENT IS POSITIVE AND THE NEGATIVE OF THE 3./5. 

POWER OF THE ABSOLUTE VALUE OF THE ARGUMENT IF THE ARGUMENT 
IS NEGATIVE. THE COMMENTS REFER TO THE POSITIVE CASE, 

THE ALGORITHM IS= 

CUBE X AND CALL THE RESULT Y, THEN WRITE Y AS 


Y = n6**C5=^N)) * (16^*M) * (Z) 

WHERE M IS BETWEEN -4 AND +4 AND Z IS BETWEEN 1/16 AND 1. 
THEN Z--**(1/5) IS APPROXIMATED BY A MINI-MAX LINEAR FIT 
FROM TWO TABLES WITH A MAXIMUM RELATIVE ERROR IN THE 
APPROXIMATION OF 5.1E-4. THEN THE INITIAL ESTIMATE OF 

T**1/5 = (16**(M/5n * (Z**1/5) 

IS REFINED BY TWO APPLICATIONS OF NEWTON’S METHOD. 

X**3/5 IS THEN CALCULATED FROM,(16**N) * (T**1/5). 


USING 

15 

B 

FIRST 

DC 

X '03' 

DC 

CL3'F35' 

DS 

18F 

STM 

14, 12, 12(13) 

L 

1,0(1) 

LD 

4,0(1) 

LR 

9, 13 

STD 

4, ARG 

LA 

13, AREA 

DROP 

15 

USING 

AREA, 13 

LTDR 

4,4' 

ST 

9,4(13) 

BHM 

NONNEG 

NI 

ARG,X'7F' 

0 

9,=X'80000000 

ST 

13,8(9) 

SR 

3,3 

IC 

3, ARG 

MV I 

ARG,X'40' 

LD 

4, ARG 

L 

4,=f '64' 

SR 

3,4 

MOR 

4,4 

M 

2,=F'3' 

LA 

6, TAB1-16 

MD 

4, ARG 

SR 

5, 5 

LA 

7, TAB2-16 

STD 

4, ARG 

IC 

5, ARG 

AR 

3,5 

SR 

2,2 

SR 

3,4 

LA 

8,TAB3+16 

BHM 

NOEXTD 

L 

2,=F '-1 ' 

0 

2,=F'5' 


TELL ASSEMBLER NEXT INST ADDR IN R15 
BRANCH AROUND NAME AND SAVE AREA 
LENGTH OF NAME 
NAME 

SAVE AREA 

SAVE CALLING ROUTINE'S GPR'S 
R1 <= ADDR ARGUMENT (X) 

F4 <= X 

R9 <= ADDR OLD SAVE AREA 
ARG <= X 

R13 <= ADDR NEW SAVE AREA 
R15 NO LONGER BASE REG 
R13 NEW BASE REG 
CHECK SIGN OF X 
LINK SAVE AREAS 

IF SIGN POSITIVE - NO FIXES OR FLAGS 
TURN OFF SIGN BIT OF ARG <= IXl 
SIGN BIT R9 ON - FLAG 
LINK SAVE AREAS 
R3 <= 0 

R3 <= EXCESS 64 EXPONENET OF ARG 

ARG <= FRACTION OF ARG 

F4 <= FRACTION OF ARG 

R4 <= HEX 40 

R3 <= EXPONENT OF ARG 

F4 <= FRACTION OF ARG **2 

R3 <= EXPONENT OF ARG*^3 

R6 <= ADDR TABLE 1-16 

F4 <= FRACTION OF ARG **3 

R5 <= 0 

R7 <= ADDR TABLE 2 - 16 

ARG <= FRACTION ARG **3 

R5 <= EXCESS 64 EXPONENT OF 

FRACTION OF ARG =^*3 

R3 <= EXCESS 64 EXPONENT OF lXl**3 

R2 <= 0 

R3 <= EXPONENT OF I X I **3 
R8 <= ADDR TABLE 3+16 
IF R3 > 0 NO SIGN EXTEND 
SIGN EXTEND FOR DIVIDE 
R2 <= EXPONENT OF T 
R3 <= N (SEE COMMENTS) 
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A 

3,=F'6S' 

R3 <= EXCESS 64 EXPONENT OF 16**N 


AR 

4.2 

R4 <= EXCESS 64 EXPONENT OF T 


STC 

4^ARG 

ARG <= T 


LO 

4, ARG 

F4 <= T 


SLL 

2,2 

R2 <= DISPLACEMENT FOR TABLE 3 


SDR 

0,0 

FO <= O.DO 


L 

5, ARG 

R5 <= 1ST 4 BYTES OF T 


IC 

4.ARG+1 

LOW ORDER BYTE R4 <= FIRST BYTE OF 

* 



FRACTION OF T 


SU 

5,8 

HIGH ORDER 3 BYTES OF R5 <= BYTES 

* 



1-3 OF FRACTION OF T 


N 

4,=N 'OOOOOOFC 

' R4 <= DISPLACE FOR TABLES 1 £ 


IC 

5,ARG+4 

R5 <= BYTES 1-4 OF FRACTION OF T 


LE 

0,0(4,75 

FO <= TABLE 2 ENTRY - SLOPE 


SRL 

5,2 

LOW 3 BYTES OF R3 <= FRACTION OF 

* 

FRACTIONAL DISPLACEMENT FOR MINI-MAN LINE 


ST 

5. FRAC 

FRAC <= FRACTIONAL DISPLACEMENT FRACTION 


MV I 

FRAC,N'40' 

FRAC <= FRACTIONAL DISPLACEMENT 


ME 

D.FRAC 

FO <= FRAC!*^SLOPE 


AE 

0,0(4,65 

FO <= MINI-MAX ESTIMATE OF Z**1/5 


LTR 

9,9 

CHECK IF X NEGATIVE 


ME 

0,0(2,85 

FO <= MIN I -MAX EST OF T^--*'1/5 


BNM 

NOSIGN 

IF X POSITIVE, NO FIXES 


0 

3,=F' 128' 

SIGN OF EXPONENT OF 16^'^N MADE MINUS 

NOSIGN 

* 

LOR 

2,0 

F2 <= MINI-MAX EST OF T'-^^l/B = EST 1 

* 

BEGIN 

TWO APPLICATIONS OF NEWTON'S METHOD 

* 

WITH 

SOME GPR FIN-UPS INTERLEAVED. 


MDR 

2,2 

F2 <= EST1 2 


STC 

3.MUL 

MUL <= SIGN(X5 * IG^-'-^N 


MDR 

2,2 

F2 <= EST 1 ** 4 


L 

1,24(95 

R1 RESTORED 


LOR 

6,2 

F6 <= EST 1 ** 4 


MDR 

2,0 

F2 <= EST 1 *Tc 5 


L 

2,28(95 

R2 RESTORED 


SDR 

2,4 

F2 <= EST 1 5 - T 


DER 

2,6 

F2 <= (EST 1 5 - TJ/EST 1 4 


L 

?, 32(95 

R3 RESTORED 


ME 

2, ONES 

F2 <= (EST 1 ** 5 - T5/5'<'EST 1 ** 4 


L 

4,36(95 

R4 RESTORED 


SDR 

0,2 

FO <= EST 2 


LOR 

2,0 

F2 <= EST 2 


MDR 

2,2 

F2 <= EST 2 ** 2 


L 

5,40(95 

R5 RESTORED 


MDR 

2,2 

F2 <= EST 2 4 


L 

6,44(95 

R6 RESTORED 


LDR 

6,2 

F6 <= EST 2 t;* 4 


MDR 

2,0 

F2 <= EST 2 5 


L 

7,48(95 

R7 RESTORED 


SDR 

2,4 

F2 <= EST 2 ** 5 - T 


L 

8,52(95 

R8 RESTORED 


DDR 

2,6 

F2 <= (EST 2 ** 5 - T5/EST 2 ** 4 


SR 

15,15 

R15 <= 0 (RETURN C00E5 


MD 

2.0HE5 

F2 <= (EST 2 ** 5 - T)/5*EST 2 ** 4 


MV I 

12(95. N'FF' 

INDICATE CONTROL RETURNED 


SDR 

0,2 

FO <= EST 3 (FINAL ESTIMATE OF T»'’*1/55 


L 

9,56(95 

R9 RESTORED 


MD 

0,MUL 

FO <= ESTIMATE OF 1 X 1 ^'*3/5 


L 

13,4(135 

R13 RESTORED 


BR 

14 

RETURN 


CNOP 

4.8 

FORCE CORRECT ALIGNMENT 

FRAC 

DC 

F 'O' 


ARG 

DC 

D '0. ' 


MUL 

DC 

N'0010000000000000' 

ONES 

DC 

N'4033333333333333' 


LTORG 



TAB1 

DC 

N '40931BB3' 



DC 

N'4099CBC3' 



so 




DC 

X 

M09F7DF5' 


DC 

X 

MOA-iygcE' 


DC 

X 

M0A8EBB4' 


DC 

X 

M0ACF139' 


DC 

X 

' 40B09F20 ' 


DC 

X 



DC 

X 

'405/2005' 


DC 

X 

'40BA214F ' 


DC 

X 

' 40BCE873' 


DC 

X 

'40BFS816' 


DC 

X 

'40C204D6' 


DC 

X 

'40C46286' 


DC 

X 

'40C6A45E' 


DC 

X 

'40CSCD17' 


DC 

X 

'40CADF05' 


DC 

X 

' 40CCDC2C ' 


DC 

X 

'40CEC64A' 


DC 

X 

'40D09EE7 ' 


DC 

X 

'40D2675C 


DC 

X 

' 40D420D>' 


DC 

X 

'40D5CC70' 


DC 

X 

'40D76B12' 


DC 

X 

'40DSF093' 


DC 

X 

' 40DAS4C8 ' 


DC 

X 

'40DC0153' 


DC 

X 

'40DD73DA' 


DC 

X 

'40DEDCF1 ' 


DC 

X 

' 40E03D ID ' 


DC 

X 

'40E194DA' 


DC 

X 

' 40E2E499 ' 


DC 

X 

' 40E42CC 1 ' 


DC 

X 

'40E56DB3' 


DC 

X 

'40E6A7C8' 


DC 

V 

A 

'40E7DB51 ' 


DC 

X 

' 40E9089C ' 


DC 

X 

' 40EA2FF0 ' 


DC 

X 

' 40EB51 SF ' 


DC 

X 

' 40EC6DB7 ' 


DC 

V 

' 40ED84A1 ' 


DC 

X 

'40EE9686' 


DC 

V 

'40EFA397' 


DC 

X 

'40F0AC04' 


DC 

X 

'4CF1AFFB ' 


DC 

X 

'40F2AFA6' 


DC 

X 

' 40F3AB2D ' 


DC 

X 

'40F4A2B6' 


DC 

X 

'40F59664' 


DC 

X 

'40F6865B' 


DC 

X 

'40F772B9' 


DC 

X 

'40F85B9D ' 


DC 

f\ 

' 40F94124 ' 


DC 

X 

'40FA236A' 


DC 

X 

'40FB02S8' 


DC 

X 

'40FBDE98 ' 


DC 

X 

' 40FCB7B1 ' 


DC 

X 

' 40FDSDEA ' 


DC 

X 

'40FE6159' 


DC 

X 

' 40FF321 1 ' 

TAB2 

DC 

X 

•3F6B5E5D ' 


DC 

X 

' 3F5B5AA7 ' 


DC 

X 

'3F4FE171 ' 


DC 

X 

' 3F4736EF ' 


DC 

X 

'3F4009F9' 


DC 

X 

' 3F3AEB91 ' 


DC 

V 

r> 

'3F366146' 


DC 

X 

' 3F328EEF ' 


DC 

X 

'3F2F4AE5' 


DC 

X 

'3F2C775S ' 


DC 

X 

' 3F29FE5F ' 
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X'3F27CF79' 

X'3F25DDFD' 

ORIGINAL PAGE IS 

DC 

DC 

jQB QUMXSXi 

DC 

X!3F242004' 


DC 

X'3F22SDB5' 


DC 

X'3F2120C2' 


DC 

X '3F1FD408' 


DC 

X'3F1EA350' 


DC 

X'3F1D8B17' 


DC 

X'3F1C8870' 


DC 

X'3F1B98DF' 


DC 

X 'SFIABA^IB' 


DC 

X'3F19EAE6' 


DC 

X'3F192921 ' 


DC 

X'3F1S73A6' 


DC 

X'3F17C948' 


DC 

X' 3F172902' 


DC 

X '3F1691EC ' 


DC 

X'3F16033Q' 


DC 

X'3F157C3C 


DC 

X ' 3FHFCTD ' 


DC 

X'3F1482E0' 


DC 

X' 3FH0F75' 


DC 

X‘3F13A19S' 


DC 

X '3F1338E1 ' 


■ DC 

X'3F12D^F‘1' 


DC 

X'3F12757A' 


DC 

X'3F121A27' 


DC 

X '3F11C2B5' 


DC 

X'3F116EE4' 


DC 

X'3F111E78' 


DC 

X'3F10D13D' 


DC 

X '3F108701 ' 


DC 

X'3F103F97‘ 


DC 

X'3EFFAD4E' 


DC 

X '3EFBS9-n ' 


DC 

X '3EF78B10' 


DC 

X ' 3EF3B0AD ' 


DC 

X ' 3EEFFS2D ' 


DC 

X ' 3EEC5FC9 ' 


DC 

X ' 3EESE5DA ' 


DC 

X '3E'E588D6' 


DC 

X ' 3EE2474D ' 


DC 

X ' 3EDF 1FE8 ' 


DC 

X'3EDC1167' 


DC 

X '3ED91A9C 


DC 

X ' 3ED63A6E ' 


DC 

X ' 3ED36FD6 ' 


DC 

X ' 3ED0B9DC ' 


DC 

X'3ECE1796' 


TABS DC 

X ' 401BDBSC ' 


DC 

XM03050C0' 


DC 

X'405472D1' 


DC 

X M0930S8C' 


DC 

X '41 100D00' 


DC 

X'4nBDBSC' 


DC 

X'4130S0C0' 


DC 

X ' 415472D1 ' 


DC 

X'419308SC' 


END 
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TESTP 


TESTP CSECT 
* 

* ROUGHLY EQUIVALENT TO THE FORTRAN CODE BELOW EXCEPT THE 

* THE FUNCTION CHINV PASSED IN THE ARGUMENT LIST 

* IS IMPLEMENTED IN LINE IN THE ASSEMBLY LANGUAGE VERSION. 

* A CALL TO THE ASSEMBLY LANGUAGE VERSION SHOULD 

* PASS A SEMI-LOGARITHMIC INTERPOLATION TABLE (CHINV(884)) 

* RATHER THAN A FUNCTION NAME. ALSO THE PARAMETER ADDRESSES 

* ARE OBTAINED FROM LOCAL STORAGE IN THE CALL TO TESTP 

* NOT FROM THE PARAMETER LIST. THE PARAMETER ADDRESSES 

* ARE INITIALIZED BY THE ENTRY POINT TESTPI IN THE ASSEMBLY 

* LANGUAGE VERSION. 

* 

* NOTE THAT THIS MEANS TESTPI MUST BE CALLED BEFORE 

* THE FIRST CALL TO TESTP OR UNPREDICTABLE ABENDS 

* WILL RESULT. 

* 

* TESTP(T,TE, J,OSIG,N,TUP,DELT,FRAC,NTAB,CHINV) 

* IMPLICIT REALMS (A-H,0-Z) 

* REAL*8 T(NTAB) ,TE(NTAB), J(NTAB),OSIG(NTAB),N(NTAB) .TUP(NTAB) 

* REALM’S C/2.997925D10/FRAC,FRACI,DELT,DELTS 

* GOTO 10 

* ENTRY TESTPI (T,TE,J,OSIG,N, TUP, DELT, FRAC,NTAB, CHINV) 

* FRACI=2.D0/FRAC 

* RETURN 

* 10 DELTS = DELT=*^FRACI 

* DO 20 1=1 ,NTAB 

* TUP(n=C*J(I)''^J(I)*OSIG(I)*N(I) 

* IF(DELTS*TUP(I).LE.TE(I))GOTO 20 

* DELTS=TE(I )/TUPCI ) 

* 20 CONTINUE 

'■« DELT=FRAC*DELTS 

* DO 20 I=1,NTAB 

* TE(n=TE(I) + TUP(I)*DELT 

* 20 T(n=CHINV(TE( I) ) 

RETURN 

* END 

* 

USING *, 15 

B TFIRST BRANCH AROUND NAME, SAVE AREA ETC. 

DC X'05' 

DC CLB'TESTP' 

ENTRY TESTPI 

USING *, 15 

B IFIRST BRANCH AROUND NAME, SAVE AREA ETC. 

DC X'06' 

DC CL7'TESTPI ' 

DS 18F SAVE AREA 

DC X'00004^10' 

DC X'00000230' 

DS 10F REGISTER STORAGE 

CHOP 0,8 FORCE DOUBLE WORD ALIGNMENT 

DC D ' 2. 997925E10' 

DC X '4519AEF8FF9F9C62' 

DC 0 ' 0. ' 

DC D '0. ' 

DC D ' 0. ' 

DC D ' 0. ' 

DC D ' 0. ' 

DC X '4000000000000000' 

FIRST ENTRY POINT 

STM 14,12,12(13) SAVE CALLING ROUTINE'S GPR'S 

LR 2,13 R2 <= ADDR OLD SAVE AREA 


TESTPI 


AREA 

TDISP 

TEND 

REGS 

C 

TION 

FRAC 

FRACI 

DELT 

WM1 

UP1 

FLOAT 

* 

* 

* 

IFIRST 
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LA 

13, AREA 

R13 <= ADDR NEU SAVE AREA 


DROP 

15 

R15 NO LONGER BASE REG 


USING 

AREA, 13 

R13 NEW BASE REGISTER 


ST 

2,4(131 

LINK SAVE AREAS 


ST 

13,8(21 



LM 

3,12.0(11 

R3-R12 <= ADDR'S ARCS 


LD 

0,0(101 

FO <= FRAC 


L 

11,0(111 

R11 <= NTAB 


LD 

2,=D'2. ' 

F2 <= 2. DO 


SLA 

11.3 

R11 <= NTAB*8 


STD 

O.FRAC 

FRAC(L0CAL1 <= FRAC 


SR 

3,4 

R3 <= BASE T - BASE TE 


DOR 

2,0 

F2 <= 2.D0/FRAC 


SR 

5,4 

R5 <= BASE J - BASE TE 


STD 

2,FRACI 

FRACI <= 2.00/FRAC 


SR , 

6,4 

R6 <= BASE OSIG - BASE TE 


SR 

7,4 

R7 <= BASE N - BASE TE 


SR ^ 

8,4 

R8 <= BASE TUP - BASE TE 


LA 

10,8 

RIO <= 8 (INCREMENT) 


SR ' 

11,10 

R11 <= 8*(NTAB-1) 


AR 

11,4 

R11 <= BASE TE + 8*(NTAB-1) COMPARAND 


STM 

3, 12, REGS 

REGS <= R3-R12 


L 

13.4(131 

R13 <= ADDR OLD SAVE AREA 


LM 

14, 12, 12(131 

GPR'S RESTORED 


SR 

15,15 

R15 <= 0 (RETURN CODE) 


MV I 

12(131, X'FF' 

INDICATE CONTROL RETURNED 


BR 

14 

RETURN 


DROP 

13 



USING 

TEST?, 15 


TFIRST 

STM 

14,12,12(13) 

SAVE CALLING ROUTINE'S GPR'S 


LR 

2,13 

R2 <= ADDR OLD SAVE AREA 


LA 

13, AREA 

R13 <= ADDR NEW SAVE AREA 


DROP 

15 

R15 NO LONGER BASE REG 


USING 

AREA, 13 

R13 NEW BASE REGISTER 


ST 

2,4(131 

LINK SAVE AREAS 


ST 

13,8(21 



LM 

3, 12, REGS 

SET UP GPR'S 


LD 

6,0(91 

F6 <= BELT 


MD 

6,FRACI 

F6 <= OELT*FRACI = DELTS 

L00P1 

LD 

0,0(5,41 

FO (= J(I) 


MDR 

0,0' 

FO <= J(I)*J(I1 


MD 

0, C 

FO <= C*J(I1*J(I1 


MD 

0,0(6,41 

FO <= C*J(I 1*J(1)-'*'0SIG(I) 


MO 

0.0(7,41 

FO <= C*J(I1>'*=J(I)^0SIG(I1*N(I1 


STD 

0, 0(8, 41 

TUP(I) <= C*J(I1^'J(I)*0SIG(I)»-'N(I) 


MDR 

0,6 

FO <= TUPd I'^-DELTS 


CD 

0,0(41 

IF(DELTSnUP(I) .LE.TE(H) 


BNH 

ARND 

GOTO ARND 


LD 

6,0(41 

F6 <= TE(I) 


OD 

6,0(8,41 

F6 <= TE(I)/TUP(I1 =- DELTS 

ARND 

BXLE 

4, 10,LOOP1 

I <= 1+1 £ GOTO L00P1 IF NOT DONE 


MD 

6, FRAC 

F6 <= DELTS''^FRAC = DELT 


L 

4,REGS+4 

R4 <= BASE TE 


STD 

6, 0^91 

DELT <= OELTS^FRAC 

L00P2 

LD 

0,0(8,41 

FO TUP(I) 


MO 

0,0(91 

FO <= TUPdl^'DELT 


AD 

0,0(41 

FO <= TEd) + TUPd)*DELT 


STD 

0,0(41 

TE(I) <- TE(I) + TUPd )-'*'DELT 


I.H 

2,0(41 

R2 <= HIGH ORDER BYTES OF TE(I) 


five 

FL0AT+1(61,2 

(41 FLOAT <= FRACTIONAL UiSPLACEMENT 


S 

2, TOISP 

REDUCE R2 BY TDISP NOW « DOUBLE WORDS 

* 



FROM BASE OF INTERPOLATION TABLES 


LD 

4, FLOAT 

F4 <= FRACTION 0 LE FRAC LT 1 


BM 

LOUT 

IF RESULT NEGATIVE - OUT OF RANGE 

* 



GOTO LOUT 


f> 

V 

2, TBNO 

COMPARE R2 TO TBND IF GREATER 


BH 

HITE 

OUT OF RANGE -- GOTO HITE 


SD 

4,=0' .5' 

F4 <= X = FRAC - .5 -.5 LE X LE .5 


SLA 

2,3 

R2 <= R2*8 NOW BYTE DISPLACEMENT 
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FROM BASE OF INTERPOLATION TABLE. 

* 

NOW 

COMPUTE WEIGHTS 

FOR CUBIC INTERPOALTION OF 


FUNCTIONS' OF T 



LDR 

2,4 

F2 <= X 


MDR 

4,4 

F4 <= X**2 = X2 


HDR 

4,4 

F4 <= X2/2 


SD 

4,=D' 1. 125' 

F4 <= X2/2 - 9/8 


LDR 

6,4 

F6 <= X2/2 - 9/8 


HDR 

4,4 

F4 <= X2/4 - 9/16 


MDR 

6,2 

F6 <= X3/2 - 9X/8 


LCDR 

0,4 

FO <= -X2/4 + 9/16 


ADR 

0,6 

FO <= X3/2 - X2/4 - 9X/8 + 9/16 


STD 

0,WM1 

WMl <= WEIGHT FOR TABLE ENTRY CORRESPOND 

* 



ING TO CLOSEST SMALLER VALUE OF T. 


LCDR 

0,4 

FO <= -X2/4 + 9/16 


SDR 

0, 6 

FO <= -X3/2 - X2/4 + 9X/8 + 9/16 


STD 

0,WP1 

WP1 <= WEIGHT FOR TABLE ENTRY CORRESPOND 

* 



ING TO CLOSEST LARGER VALUE OF T. 


ADR 

6,2 

F6 <= X3/2 - X/8 


AD 

4,=D'.5' 

F4 <= X2/4 - 1/16 


MD 

6,=X'4055555555555555' F6 <= X3/6 - X/24 : 


LDR 

0,4 

FO <= X2/4 - 1/16 1 


SDR 

0,6 

FO <= -X3/6 + X2/4 + X/24 - 1/16 


ADR 

6, 4 

F6 <= X3/6 + X2/4 - X/24 - 1/16 


MD 

0,0(12,2) 

FO <= WEIGHT1 * CHNV1 


MD 

6,24(12,2) 

F6 <= WEIGHT4 * CHNV4 


LD 

2,WM1 

F2 <= WWIGHT2 


ADR 

0,6 

FO <= W1*CHNV1+W4*CHNV4 


LD 

4,UP1 

F4 <= WEIGHTS 


MD 

2,8(12,2) 

F2 <= W2*CHNV2 


MD 

4, 16(12,2) 

F4 <= W3*CHNV3 


ADR 

0,2 

FO <= W1*CHNV1+W2*CHNV2+W4^CHNV4 


ADR 

0,4 

FO <= INTERPOLATED VALUE OF CHINV(TE(D) 


STD 

0,0(3, 4) 

T(I) <= GHINVCTEd)) 


BXLE 

4, 10, L00P2 

K=I + 1 £ GOTO L00P2 IF NOT DONE 


B 

ARND1 


LOUT 

STD 

0,0(3, 4) 

T(I) <= TE(I) (FOR T < 4096 K) 


BXLE 

4, 10, L00P2 

I <= 1+1 AND GOTO L00P2 IF NOT DONE 


B 

ARND1 

GOTO ARND1 

HITE 

SD 

0,TI0N 

FO <= TE(I) - TION 


HDR 

0,0 

FO <= (TE(I)-TI0N)/2.D0 


STD 

0,0C3,4)- 

T(I) <= (TE(I )-TI0N)/2. DO 


BXLE 

4, 1Q,L00P2 

I <= 1+1 AND GOTO L00P2 IF NOT DONE 

ARNDT 

L 

13,4(13) 

R13 <= ADDR OLD SAVE AREA 


LM 

14,12,12(13) 

GPR'S RESTORED 


SR 

15, 15 

R15 <= 0 (RETURN CODE) 


MV I 

12(13), X'FF' 

INDICATE CONTROL RETURNED 


BR 

14 

RETURN 


END 
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DRIGINSE PAGE IS 
:0E P(X)R QUAUTt 

TOUT 


lioUT CSECT 

^ ROUGHLY EQUIVALENT TO THE THREE FORTRAN SUBROUTINES BELOW 

^ NOUT OPENS LOGICAL UNIT 9 CFT09F001) AND DOES THE OUTPUT 

i OF THE SUBROUTINE NOUT. THE PARAMETER NTAB IS PASSED 

^ BY ENTRY POINT NOUT AND THE NTAB IN THE CALLING SEQUENCE 

^ TO TOUT IS IGNORED BY THE ASSEMBLY LANGUAGE VERSION OF 

THESE ROUTINES. FORTRAN CLOSES DATA SETS IT KNOWS ABOUT 
^ • BUT FORTRAN WON'T KNOW ABOUT THIS DATA SET SO CTOUT MUST 

iF BE CALLED BEFORE THE STOP STATEMENT IN THE MAIN ROUTINE. 

^ THE FORTRAN CTOUT DOES THE SAME THING - I.E. IT CAUSES A 

4: CLOSE TO BE ISSUED FOR THE DATA SET REFERENCED BY THE 

t DDNAME FT09F001. USES QSAM UNDER 0S/VS2. 

* 

NOTE S:^*'*-'^^*^'*^''*:*'*:**^*;*-**'*:;*;^^'^’****'*:*;*:*** 


:L'F** NOUT must BE CALLED BEFORE TOUT IF TOUT IS TO 

WORK SINCE NOUT SETS UP AN AREA WITH THE 
*:«.■** REGISTERS THAT TOUT USES - THIS MAKES SENSE 

IN this APPLICATION (NOUT IS ALWAYS CALLED 
FIRST) BUT MUST BE CHANGED IF NOUT IS NOT TO BE 
**** CALLED BEFORE TOUT 


***;*;*:**s:*:*:******.-*:****:*:*.'*:*»:**:*:**;»:**3;********^-^^^'^“K*.1;***»:***.i:*^^-:*;.-i: 

* 

* 

SUBROUTINE NOUT ( EF LUX , PS I 0, FRAC , T I MMAX , NTAB , N , S ) 

IMPLICIT REAL-^^S (A-H,0-Z) 

REALMS N(NTAB) ,S(NTAB) 

FORMAT (10A8) 

WRITEC9, 9001) EFLUX, PSIO,FRAC,T I MMAX, NTAB 
WRITE(9,°001)N 
WRITE(9,. 9001 )S 
RETURN 
END 

SUBROUTINE TOUT (TI M, T , J, D ELT, II TER , NTAB ) 

IMPLICIT REAL*fi (A-H,0-Z) 

REAL=*-'8 T(NTAB) , JCNTAB) 

FORMAT! 10A8) 

WRITE(9,9001)TIM,DELT,IITER 
WRITE(9,9001)T 
WRITE(9,9001)J 
RETURN 
E^ND 

SUBROUTINE CTOUT 
END FILE 9 
RETURN 
END 

USING *,15 
B WFIRST 

DC X'04' 

DC CLS'TOUT ' 

ENTRY NOUT 
USING *, 15 

NOUT B OFIRST 
* 


BRANCH AROUND NAME, VARIABLES, SAVE AREA 
AND OTHER ENTRY POINTS 
LENGTH OF NAME 
NAME 


BRANCH AROUND NAME, VARIABLES, SAVE AREA 
AND OTHER ENTRY POINT 


* 

* 

*9001 

*: 

*• 

* 

* 

* 

* 

* 

*9001 

* 

* 

* 

.* 

* 

% 


86 




DC 

X '04' 

LENGTH OF NAME 


DC 

CLS'NOUT ' 

NAME 


USING 

*, 15 



ENTRY 

CTOUT 


CTOUT 

B 

CFIRST 

BRANCH AROUND NAME, VARIABLES AND AREA 


DC 

X '05' 

LENGTH OF NAME 


DC 

CLS'CTOUT ' 



DROP 

15 



USING 

AREA, 13 


AREA 

DS 

18F 

SAVE AREA 

BUF 

DS 

20F 

OUTPUT BUFFER (ONE CARD - WE USE QSAM) 

BLNCRD 

DC 

20CL4' ' 

A BLANK CARD - I'M LAZY 

TEN 

DC 

F ' 10' 


LM0V3 

MVC 

BUF(1),0(3) 

TO BE EXECUTED BY AN EX 

LMOV-4 

nvc 

BUF(1),0(4) 


SREG 

DS 

5F 



DROP 

13 



USING 

CTOUT, 15 



CTOUT 



CFIRST 

STM 

14,12,12(13) 

SAVE CALLING ROUTINE'S GPR'S 


LR 

2, 13 

R2 <= ADDR OLD SAVE AREA 


LA 

13, AREA 

R13 <= ADDR NEW SAVE AREA 


DROP 

15 

R15 NO LONGER BASE REG 


USING 

AREA, 13 

R13 NEW BASE REG 


ST 

2,4(13) 

LINK SAVE AREAS 


ST 

13,8(2) 



CLOSE 

LU9DGB 



RETURN 

1 SEQUENCE 



L 

13,4(13) 

R13 <= ADDR OLD SAVE AREA 


LH 

14,2, 12(13) 

GPR'S RESTORED 


SR 

15, 15 

R15 <=0, RETURN CODE 


MV I 

12(13), X'FF' 

INDICATE CONTROL RETURNED 


BR 

14 

RETURN 


DROP 

13 


* 

NOUT(EFLUX,PSIO,FRAC 

;,TIMMAX,NTAB,N,S) 


USING 

NOUT, 15 


OFIRST 

STM 

14, 12,12(13) 

SAVE CALLING ROUTINE'S GPR'S 


LR 

2, 13 

R2 <= ADDR OLD SAVE AREA 


LA 

13, AREA 

R13 <= ADDR NEW SAVE AREA 


DROP 

15 

R15 NO LONGER BASE REG 


USING 

AREA, 13 

R13 NEW BASE REG 


ST 

2,4(13) 

LINK SAVE AREAS 


ST 

13,8(2) 



LM 

2, 8, 0(1) 

R2-R8 <= ADDR 'S ARGS 


OPEN 

(LU9DCB, (OUTPUT) ) 


MVC 

BUF(80), BLNCRD 

1 ' FILL BUFFER WITH BLANKS 


MVC 

BUF(8),0(2) 

FIRST 8 BYTES OF BUFFER <= EFLUX 


MVC 

BUF+8(8),0(3) 

2ND 8 BYTES OF BUFFER <= PISO 


MVC 

BUF+16(8),0(4) 

3RD 8 BYTES OF BUFFER <= FRAC 


MVC 

BUF+24(8) , 0(5) 

4TH I BYTES OF BUFFER <= TIMMAX 


MVC 

BUF+36(4),0(6) 

2ND HALF OF 5TH 8 BYTES OF 




BUFFER <= NTAB 


PUT 

LU9DCB, BUF 

WRITE OUT 1ST RECORD 


LR 

3,7 

R3 <= BASE ADDR N 


LR 

4,8 

R4 <= BASE ADDR S 


L 

7, 0(6) 

R7 <= NTAB 


SR 

6 , 6 

R6 <= 0 


D 

6, TEN 

R6 <= REMAINDER OF NTAB/10 




R7 <= INTEGER PART OF NTAB/10 


LR 

5, 6 

R5 <= REMAINDER OF NTAB/10 


M 

6, TEN 

R7 <= (NTAB/10)*10 INTEGER MODE 




# CARDS - 1 (UNLESS NTAB ENDS IN 0) 


SLL 

5,3 

F5 <= REMAINDER NTAB/10 * 8 
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SLL 

7,3 

R7 <= (NTAB/10)*10 * 8 

* 



« BYTES IN # CARDS - 1 


LA 

6,80 

R6 <= 80 (INCREMENT FOR LOOP) 


LR 

2,3 

R2 <= BASE ADDR N 


LA 

4,0(4) 

HIGH BYTE OF BASE ADDR S ZEROED 


SR 

7.6 

R7 <= (NTAB/10 - 1) 80 


STM 

5, 7.SREG 

SAVE INCREMENTS AND OFFSET FOR COMPARAND 


STM 

6,7,SREG+12 

TWO COPIES OF INCREMENTS AND OFFSET 


AR 

7,3 

R7 <= BASE ADDR N + ( (NTAB/1 0)- 1 )=^S0 

L00P1 

MVC 

BUF(SO) ,0(3) 

PUT NEXT 80 BYTES IN BUFFER 




AND WRITE THEM OUT 


PUT 

LU9DCB.BUF 



BXLE 

3,6,L00P1 

KEEP GOING TILL WE'RE DONE 


LTR 

5,5 

CHECK IF NO MORE THINGS TO WRITE 


B2 

ARNDl 

IF NOT DON'T WRITE OUT ANOTHER CARD 

* 



BECAUSE FORTRAN WOULDN'T 


MVC 

BUF(80),BLNCRD 

1 FILL BUFFER WITH BLANKS 


EX 

5, LM0V3 

MOVE IN LAST FEW (<10) VALUES 

* 



AND WRITE OUT LAST CARD FOR N 


PUT 

LU90CB,BUF 


ARND1 

SR 

7,2 

R7 (= (NTAB/10 - 1) * SO 


AR 

7,4 

R7 <= BASE S + ( (NTAB/10)-1 )*80 

L00P2 

MVC 

BUF(80),0(4) 

PUT NEXT 80 BYTES IN BUFFER 




AND WRITE THEM OUT 


PUT 

LU9DCB, BUF 



BXLE 

4,6,L00P2 

KEEP GOING TILL WE'RE DONE 


LTR 

5,5 

CHECK IF NO MORE THINGS TO WRITE 


B2 

ARND2 

IF NOT DON'T WRITE OUT ANOTHER CARD 

* 



BECAUSE FORTRAN WOULDN'T 


MVC 

SUF(80),BLNCRD 

I FILL BUFFER WITH BLANKS 


EX 

5,LM0V4 

MOVE IN LAST FEW (<10) VALUES 

* 



AND WRITE OUT LAST CARD FOR S 


PUT 

LU9DCB, BUF 


* 

* 

RETURN 

1 SEQUENCE 


ARND2 

L 

13,4( 13) 

R13 <= ADDR OLD SAVE AREA 


LM 

14,8, 12(13) 

GPR'S RESTORED 


SR 

15,15 

R15 <= 0, RETURN CODE 


MV I 

12C13),X'FF' 

INDICATE CONTROL RETURNED 


BR 

14 

RETURN 


DROP 

13 



TOUT(TIM,T,J,DELT/HTER,NTAB) WE IGNORE LAST PARAMETER 


USING 

TOUT, 15 


UFKRST 

STM 

14,12,12(13) 

SAVE CALLING ROUTINE'S GPR'S 


LR 

2,13 

R2 <= ADDR OLD SAVE AREA 


LA 

13, AREA 

Rl3 <= ADDR NEW SAVE AREA 


DROP 

15 

R15 NO LONGER BASE REG 


USING 

AREA, 13 

R13 NEW BASE REG 


ST 

2,4(13) 

LINK SAVE AREAS 


ST 

13,8(2) 



LM 

2.6,0( 1) 

R2-R6 <= ADDR 'S ARG'S WE USE 


MVC 

BUF(80) , BLNCRD 

1 . FILL BUFFER WITH BLANKS 


MVC 

BUF(8) ,0(2) 

FIRST 8 BYTES OF BUFFER <= TIM 


MVC 

BUF+8(8) , 0(5) 

2ND 8 BYTES <- DELT 


MVC 

BUF+20(4),0(6) 

2ND HALF 3 8 BYTES <C= IITER 


PUT 

LU9DCB, BUF 



LM 

5, 9,SREG 

R5-R9 <= NUMBER OF VALUES ON LAST CARD FOR 

* 



T AND J AND INCREMENTS AND COMPARANDS 


AR 

7,3 

R7 <= COMPARAND FOR L00P3 


AR 

9, 4 

R9 <= COMPARAND FOR LOOP4 

L00P3 

MVC 

BUF(80) ,0(3) 

PUT NEXT SO BYTES IN BUFFER 

* 



AND WRITE THEM OUT 


PUT 

LU9DCB, BUF 



BXLE 

3, 6, L00P3 

KEEP GOING TILL WE'RE DONE 


LTR 

5,5 

CHECK IF NO MORE THINGS TO WRITE 


B2 

L00P4 

IF NOT DON'T WRITE OUT ANOTHER CARO 
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* 



BECAUSE FORTRAN WOULDN'T 



MVC 

BUF(SO) , BLNCRD FILL BUFFER WITH BLANKS 


EX 

5, LM0V3 

MOVE IN LAST FEW (<10) VALUES 





AND WRITE OUT LAST CARD FOR T 



PUT 

LU9DCB, BUF 



L00P4 

MVC 

BUF(80),0(4) 

PUT NEXT 80 BYTES IN BUFFER 


* 



AND WRITE THEM OUT 



PUT 

LU9DCB,BUF 




BXLE 

4, 8, LOOP4 

KEEP GOING TILL WE'RE DONE 



MVC 

BUF(SO) , BLNCRD FILL BUFFER WITH BLANKS 


LTR 

5,5 

CHECK IF NO MORE THINGS TO WRIT 

E 


BZ 

ARND3 

IF NOT DON'T WRITE OUT ANOTHER 

CARD 

* 



BECAUSE FORTRAN WOULDN'T 



EX 

5,LMOV4 

MOVE IN LAST FEW (<10) VALUES 


* 



AND WRITE OUT LAST CARD FOR J 


* 

PUT 

LU9DCB, BUF 



* 

* 

RETURN 

SEQUENCE 



ARND3 

L 

13,4(13) 

R13 <= ADOR OLD SAVE AREA 



LM 

14,9, 12(13) 

GPR'S RESTORED 



SR 

15,15 

R15 <= 0, RETURN CODE 



MV I 

12( 13),X 'FF' 

INDICATE CONTROL RETURNED 



BR 

14' 

RETURN 


LU9DC8 

DCB DEVD 

=DA,MACRF=PM, 

DSORG=PS,RECFM=FB, LRECL=80, DDNAME=FT09F00 


ENO 






Appendix B 


STEADY STATE MODEL OF THE SOLAR ATMOSPHERE 

We have constructed a steady state numerical model of the solar 
atmosphere. The model was developed to investigate the effects of upward 
velocities and diverging magnetic field patterns on the temperature and 
density structure of the solar atmosphere; however, for this work the 
model is used only to provide reasonable temperature and density profiles 
for the estimation of the effect of reverse current heating on the atmo- 
sphere. The computer program calculates the run of temperature and 
density in an individual flux tube. 

The equations governing the behavior of an inviscid compressible 
fluid in the presence of gravity are 

1^ + V • (p-3) = 0 , (B.i) 

(pu) + V • (p'Sit) = - VP + pg , (b.2) 

^ — * — » — + 

(pe) + V • (pcu) =-V*q-P7*u-£+S, (b.3) 

where e is the total internal energy per unit mass, q is the heat 
flux, g is the gravitational acceleration, £ is the energy lost via 
radiation, and S is the sum of all other non-thermal energy sources or 
sinks. For flow along a magnetic flux tube, considering variation only 
along the field lines and assuming the radius of curvature of the field 
lines to be large compared to the dimensions of the flux tube, we see 
that the equations become one-dimensional. If we add the definition of 
the heat flux and an equation of state to Equations (b.1)-(b. 3), we may 
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wl'ite a coiuplGtc sot of equations foi’ the steady state (3/^t"0) case 


fiP 

ds 


•1- p 


du 

£ 

ds" 




dA 

ds 


( B .!0 


du 

H 

o P — 
s ds 


pu 


d€ 




s ds ds 


- P 


du 

g 

ds 


dP 


A ds 


(baO 


,, dT 

q 1C — . 

*s ds ■’ 


P - pKT , 
u ' 


(b.t) 

(b.8) 


where A is the ai’ea of the flux tube vi is the mean particle luass^ 

1C is the heat conductivitj'', k is BoltKinann's constant^ T is the 
fluid teinperaturoj s measxiros distance along' the flux tube and the 
subscript s denotes the component of a vector along' the flux tube. We 
have neglected transport of energy and momentum acros.s field lines in 
writing oqvxations (B.)|) -(b. t) . We wish to apply Equations (b.J|-)-(b. 8) 
to the solar atmosphere. Eor this case we shall assvime the plasma to 
be pure hydrogen except for compufing the radiative losses. To account 
for radiative losses, wo have assumed that it is reasonable to treat the 
solar atmosphere as optically thin (we discuss this assumption later). 

We have adopted the radiative loss function calculated by Raymond e_t al. 
(I97w) as modified by Raymond (l9‘?h) to include radiative losses from 
Ar and neutral hydrogen excitation, but excluding radiative losses duo 
to forbidden linos for temperatures below T . We have used the 
values of u and 1C derived by Moore and Fung (Ixl'f;’?) foi' a pure hydro- 
gen plasma. With these and the oqxmtion of state, we may el'iminate the 



pressure. We choose as our dependent variables q ^ T , u and n , 

s s 

the number density of hydrogen nuclei and rewrite (b.4)-(b.y) in a form 
more convenient for numerical solution; 



(b-9) 


dn 

ds 



~ s dA dT (Itx)nk ^ nkT dX 

A ds ds m m„ dT 

n H 


r 

(b.io) 


£ + S - ,4- 'Be — 

2 s ds 


(l*x) . (T-T^) § + (l-,-X).A^ ^ 


(B.ll) 


nu^A =: constant (B.12) 

where m^^ is the mass of a hydrogen atom, X is the fraction of hydrogen 
nuclei that are ionized, and is the hj'drogeu ionization energy 

expressed as a temperature, l.Ql!’? x ^0 K. Equation (b. 12) is the inte- 
gral of Equation (b.^I-), we need only solve three first order ordinary 
diff ei’ential equations to calculate the run of temperature and density 
in a flux tube. 


We have written an assembly-language subroutine, to execute on IBM 

♦ dT dn 

360 or 370 series computers, to evaluate the quantities — , ^ — and 

dq ® 

s dA 

—— , given A , — , g , S , n , T and q , This subroutine may 
ds ds s s 

be used with a standard library ordinary differential solvei’, or as we 

have done with one coded specially for this problem. The quantities 
dX dA 

g , } ""d S are tabulated as a function of T 

Cl u S Cl S 
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and s , and the values for a particular T or s are computed by a 
cubic interpolation scheme similar to the one described in Appendix A. 

The subroutine we have written to solve the coupled set of ordinary 
differential equations (b.9)-(b. 12) uses an Adams-Bashf orth-Moulton 
fourth order lineai* multistep integration scheme (see Isaacson and Keller 
1966 ) with a fourth order Runge-Kutta scheme (with a smaller step size) 
to "start up" the linear multistep method and provide intermediate values 
when halving the step size. The routine returns the values of T , u , 
q and n at intervals from the starting point Specified by the calling 
program and reduces the step size or increases it according to the 
requested accuracy. The pretabulated quantities are read in by the main 
program which also reads in starting values, calls the differential 
equation solver and writes out the results of the integration. 

The downward heat flux in the corona above an active region is 

~ 9 ^ 10^ (Noyes 197^) • Since the thermal conductivity of the solar plasm 

*5 /o 

is a strong function of temperature (“ T^ ) , this heat flux must be 
largely radiated away above the low chromosphere. We have the choice of 
starting with our initial values where the heat flux is large (in the 
corona) and calculating the solutions to a region where the heat flux is 
small (the chromosphere), or proceeding in the reverse direction from 
the region whei'e the heat flux is small, it is well known that the 
latter choice is preferable mmierically (Acton I97O, Isaacson and Keller 
1966). This is basically because the numerical calculation proceeding 
from the region of large heat flux to the region of small heat flux is 
not a "well posed" problem ( Isaacson and Keller I966) since a small 
relative change in the initial value of the heat flux can cause a large 
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relative change in the final value. We therefore shall choose our 
starting point near the temperature minimum. 

There are three major difficulties with starting the calculation 

4 

below 3 ^ 10 K. The first is that the atmosphere becomes optically thick 
and therefore the radiative losses cannot be calculated simply. Second, 

the radiative loss function calculated by Raymond is not tabulated below 

4 

10 K. Third, the approximation that the atmosphere is purely hydrogen 
breaks down as the fraction of ionized hydrogen becomes very small 
because the electron density (which appears in the expression for the 
radiative losses) is grossly underestimated by (a. 4), since the major 
contribution to the electron density is from trace elements with low 
ionization potentials (e.g. Na) . However, for the purposes of this work, 
we only need a model that represents the overall structure of the atmo- 
sphere reasonably well. This is particularly true since (cf . Chapter 3) 
the calculation of the heating of the cool dense portions of the atmo- 
sphere by the reverse current is not accurate after the first few tenths 
of a second due to the neglect of Coulomb collisions. We do not attempt 
a solution of the radiative transfer problem. We use a power law extra- 
polation of Raymond's (I 976 ) radiative loss coefficient. We also use 
Equation (A. 4) to find the electron density. The fact that the atmosphere 
is not optically thin is compensated for by the underestimate of the 
electron density. We have extrapolated Raymond's (I976) radiative loss 
function with a power law above and below the tabulated range 

g 8 

(t=10 -T 10 k) . For the high temperatures above T=10 K , this should be 
a reasonable approximation since the losses for these temperatures are 
almost completely due to thermal bremsstrahlung and therefore should vary 
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1/p 

as ^ T ; howevex/ these tempex’atures are not o± importance in the 

ll- 

present calculation. The powei- law extrapolation below 10 K is purely 
ad hoc, but the range over which extrapolated values are used is small 
a factor of 2) and the calculation of X’adiative losses for these 
temperatures is at best approximate in any event. The resulting tempera- 
tvu’e and density profiles resemble the solar atmosphere in overall 
structure. Since the atmosphere varies from active region to active 
region, this should provide an adequate representation lor the purposes 
of the calculations of Chapter 3. 

To produce the model used (see Chapter 3)^ we integrate up from 
near the temperature minimum (T~i^200K, n=l.lQ21!' X 10^^). The heat flux 
and velocity are taken to be zero at this point. No non-thennal enei’gy 
input was included in the calculation. The resulting temperature, den- 
sity and heat flux at the top of the model (corresponding to the injection 

point for the beam in Chapter 3) were Ts3 X 10*^, nal X 10^ cm~^ and 
o -2 -1 

Fsfci .30 X 10 ei’g cm s ^ iu reasonable agreement with the values given 
by Noyes (l9Tl) • 

Listings of two main programs and several subroutines are pi^ovided 
for the sake of completeness. The first main program and associated 
subrovitines pi'odvice the tables that are i*eq\iii'ed for the cubic intei’pola- 
tion. The second main program reads in starting valiies lor the solution 
of the coupled set of differential equations and writes out the results 
both as tables suitable foi' people to look at and (if desired) for 
machines to read. The sxibroutine ABMINT is the differential equation 
solver described above. The present version is in I’ORTK.'^ and is cer- 
tainly adequate lor the purpose of this work. An adaptation of the 



present main program to solve a boundary value problem rather than an 
initial value problem would (absent the wealth of Croesus) require this 
routine to be hand coded. The assembly language subroutine DIVF calcu- 
lates the quantities needed bj'^ ABMINT to integrate the differential 


equations . 
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TABULATION ROUTINES 


IMPLICIT REALMS (A-H,0-Z) 

REALMS G(820I, DADS (-320) ,S0R( 820) ,ONEK( 8201 , CHI (820) , DCHI ( 820) , 
,A(820) , LUM ( 820) jTST/Z-I'I 10000000000000/, DT/Z43 10000000000000/, 
•SST/Z 47 1000 0000000000/, DS/Z4 6 10000000000000/ 

COMMON /PARAM/ FRAC , AMP , SCALE 
9001 FORMATdOAS) 

5001 F0RMAT(r'13.6) 

THIS PROGRAM CALCULATES SEMI-LOGARITHMIC INTERPOLATION 
TABLES FOR A STEADY STATE MODEL OF THE PLASMA IN A MAGNETIC 
FLUX TUBE UNDER THE INFLUENCE OF GRAVITY. FOUR FUNCTIONS 
OF TEMPERATURE AND FOUR FUNCTIONS OF S (DISTANCE ALONG THE 
FLUX TUBE FROM THE SUN'S SURFACE) ARE TABULATED. THE 
FUNCTIONS OF TEMPERATURE (T) ARE= 

THE INVERSE OF THE THERMAL CONDUCTIVITY (ONEK) 

THE RADIATIVE LOSS COEFFICIENT (LUN) 

THE IONIZATION FRACTION (CHI) 

THE DERIVATIVE OF THE IONIZATION FRACTION (DCHI) 

THE FUNCTIONS OF S ARE: 

THE FORCE OF GRAVITY ALONG THE TUBE (G) 

THE AREA OF THE TUBE (A) 

THE LOGARITHMIC DERIVATIVE OF THE AREA (DADS) 

THE NON-THERMAL ENERGY INPUT (SOR) 

THE TABULATION RANGE IN TEMPERATURE IS 4.096E3 - 6.71E7 (K) . 

THE TABULATION RANGE IN S IS 1.677E7 - 2.75E11 (CM). 

THE TABLES ARE WRITTEN OUT TO FORTRAN LOGICAL UNIT 9 AND THE 
PROGRAM READS IN 641 VALUES OF TEMPERATURE AND RADIATIVE 
LOSS COFFICIENT (RAYMOND, PRIVATE COMMUNICATION) USED TO 
TABULATE THE RADIATIVE LOSS COFFICIENT. 

THE PROGRAM ALSO READS IN SEVERAL PARAMETERS THAT CHARACTERIZE 
THE FLUX TUBE AND THE NON-THERMAL ENERGY INPUT: 

FRAC: THE AREA OF THE FLUX TUBE IS ((D+S)/D)**2 WHERE 

D IS FRAC TIMES A SOLAR RADIUS. 

AMP: THE INTEGRAL OF THE NON-THERMAL ENERGY DEPOSITED IN 

A FLUX TUBE OF CONSTANT AREA IS AMP (ERG PER CM^^^^Z PER SEC). 

scale: the form OF THE NON-THERMAL ENERGY INPUT IS 

(COS( (S'-FPI)/(2*SCALE) ) ):«:^2 FOR S LESS THAN SCALE 
AND ZERO FOR S GREATER THAN SCALE. SCALE IS INPUT 
IN SOLAR RADII (INPUT OF 1. MEANS SCALE IS ABOUT 7.E10 CM). 


INITIALIZE tables: 

DO 1 1=1,820 
G(I )=0. DO 
DADS( I ) =0. DO 
SORd )=0. DO 
A ( I ) = 1 . D 0 
ONEK(I)=O.DO 
CHI (I )=0.D0 
DCHI (I )=0. DO 
1 LUM(I)=O.DO 

READ IN PARAMETERS 

READ(5, 5001 )FRAC 
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READ (5,500 DAMP 
READ(5,500DSCALE 
C 

C CALCULATE FUNCTIONS OF TEMPERATURE: 

C 

DO 20 1=1,3 
T=TST-DT 
K = (I-D=k256+1 
DO 10 J=1,2-13 
CHI (K)=FCHI (T) 

DCHI(KI=FDCHI(T) 

ONEK(K)=FKAP(T) 

LUM(K)=FLUM(T) 

T=T+0T 

10 K=K+1 

TST=16. DOnST 
20 DT=16.DO^OT 
T=TST-DT 
DO 30 K=769,S2D 
CHI CK)=FCHI (T) 

DCHI (K)=FDCHI CT) 

ONEK(K)=FKAP(T) 

LUM(K)=FLUMCT) 

30 T-T+DT 
C 

C CALCULATE FUNCTIONS OF S: 

C 

DO 50 1=1 >3 
S=SST-OS 
K=(I-1)^256+1 
DO -10 J=l,2^3 
G(K)=FG(S) 

A(K)=FA(S) 

DADS(KD=FDADS(S) 

SORCK)=FSOR(S) 

S=S+DS 

^0 K=K+1 

SST=16.D0‘»-SST 
50 DS=16.D0'^DS 

S=SST'-DS 
DO 60 K=769,S20 
G(K)=FG(SI 
A(K)=FACS) 

DADS(K)=FDADS(S) 

S0R(IO=FS0R(S) 

60 S=S+OS 
C 

C WRITE OUT TABLES: 

C 

WRITECD, 9Q0DG 
WRITE(9,90QDDADS 
WRITE(9.900DSOR 
WRITE(9,900DA 
URI TE(9, 900 DONEK 
WRITE(9, 900DLUM 
WRITEC9.90QDCHI 
WRITE(9,900DDCHT 
STOP 
END 

REAL FUNCTION FCHI*8(T) 

C 

C THIS FUNCTION CALCULATES THE lONIHATlON FRACTION AS A FUNCTION 

C OF THE temperature (T). THE IONIZATION FRACTION (FCHI) IS 

C DEFINED AS NE/CHH+NPI WHERE NE IS THE NUMBER DENSITY OF 
U ELECTRONS, AND NH AND NP ARE THE NUMBER DENSITIES OF HYDROGEN 

C ATOMS AND PROTONS RESPECTIVELY, SEE MOORE AND FUNG, SOLAR 

C PHYSICS 23 ( 1973Y,7fi-|02 FOR FORMULAE. 

C 

IMPLICIT REALMS (A-H.O-Z) 
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DATA 0NE3/ZC0S5555555555555/ 

COMMON BETA, EBETA, B 1 3, TEMP 1 , TEMP2, TCHI , D 
BETA=1 . 58D5/T 
EBETA=DEXP(BETA) 

B13=SETA*Toh-E3 

TEMP1=0. 42SSD0+0. 5DO-*:DLOG(BETA) + . 46980 0*613 
TEMP2=2. 22D-6*BETA*TEMPl*EBETA 
TCHI=1 . D0/(1 . D0+TEMP2) 

FCHI=TCHI 
RETURN 
END 

REAL FUNCTION FDCHI*8(T) 

THIS FUNCTION CALCULATES THE DERIVATIVE OF THE IONIZATION FRACTION 
CD DHI / DT) AS A FUNCTION OF TEMPERATURE (T). SEE FUNCTION FCHI. 

IMPLICIT REAL*8 CA-H,0-Z) 

COMMON BETA, EBETA,B13, TEMPI, TEMPO, TCHI,D ' 

FDCHI = 1 . 406D-1 1*TCHI*TCHI*BETA*BETA*EBETA*m . D 0+BETA) *TEMP 1 
.+ ( . 500-. 1566D0*B13)D 
RETURN 
END 

REAL FUNCTION FKAP*SCT) 

THIS FUNCTION CALCULATES THE INVERSE OF THE TOTAL THERMAL 
CONDUCTIVITY AS A FUNCTION OF TEMPERATURE. THE DEPENDENCE OF THE 
CONDUCTIVITY 0N THE "COULOMB LOGARITHM" IS APPROXIMATED IN A 
MANNER SIMILAR TO MOORE AND FUNG, SOLAR PHYSICS 23 (1972), 78-102. 

IMPLICIT REALMS CA-H.O-Z) 

COMMON BETA, EBETA, B 1 3 , T EMP 1 , TEMPO , TCH I , D 
REALMS PO/0. DO/, CL/0. DO/, CK 1 /O . DO/, RKAY/ 1 . 3S062D-16/, 

.MFROT/1 .673520-24/, ESU/4. 80325D-10/,PI/Z413243F6A8885A30/ 
IFCPO.NE.O.OO)GOTO 10 
P0=1 . D5*1 . D 10*(2. D0*RKAY) 

CL=(3. D0*RKAY**2)/(DSQRT(2. D0*PI*P0)*ESU**3) 

CK1=(9. DO*RKAY*DSQRT(RKAY) )/(4. D0*DSQRT(MPR0T) ) 

10 CLAM=(T*T*CL)*OSQRT((1.DO+TCHI)/C2.00*TCHI)) 

T12=DSQRTCT) 

IFCT.GT.4.2D5)CLAM=CLAM*6.480741D2/T12 
TEMPK = (CK 1*T)/(9. 1 2D- 1 4+7 . 95D- n/CTEMPC*T 1 2) ) 

FKAP=1 .DO/CTEMPK+C 1 . S9D-5*T*T*T 1 2)/DL0G ( CLAM) ) 

RETURN 
END 

REAL FUNCTION FLUM*8CT) 

THIS FUNCTION CALCULATES THE RADIATIVE LOSS COEFFICIENT 
SUCH THAT THE RADIATIVE LOSSES FROM AN OPTICALLY THIN 
PLASMA OF SOLAR ABUNDANCESARE FLUM*(NE**2) WHERE NE IS 
THE ELECTRON NUMBER DENSITY. THE CALCULATION OF THE 
RADIATIVE LOSS COEFFICIENT IS RAYMOND'S (PRIVATE COMM.) 

IMPROVEMENT OF THE CALCULATIONS OF RAYMOND, COX AND SMITH 
AP, J. 204 (1976), 290-292. 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 T,L,TD(641),LD(641),LOGT,ERR,FINT(10),XDIF(10),WRK(10) 
REAL*4 RT0(641),RL0(641) 

LOGICAL SORT/. FALSE./, EXTRAP/ . FALSE ./, FI RST/. TRU E . / 

EQUIVALENCE (TD ( 32 1 ) > RTO ( 1 ) ) , ( LD ( 32 1 ) , RLD ( 1 ) ) 

8001, FORMAT(20A4) 

IF(FIRST)GOTO 100 
1 10 LDGT=DL0G10(T) 

IF(L0GT.LT.TD(1))G0T0 200 
IF(LOGT.GT.TD(NRADPT))GOTO 300 
ERR=-1 .DO 

CALL AITKEN(L, LOGT, 10, ERR, TD, LD, NRABPT, SORT, EXTRAP, F I NT, XD I F , WRK , 
.810,8400,8400) 

10 FLUM=10. D0**L 
RETURN 
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100 NRADPT=641 

CALCULATING LEAST SQUARE FITS FOR POWER LAW EXTENSION 
OF CALCULATED RADIATIVE LOSS COEFFICIENT BEYOND TABULATED 
RANGE. ONLY DO ON FIRST CALL. 


READ(S,8001)RTD 
READ(S,8001)RLD 
DO 15 I=1,NRADPT 
TD( n=OBLE(RTD(I) ) 
15 LD(n=OBLE(RLD(I ) ) 
A1=0. DO 
B1=0.D0 
TEMP1=O.DO 
TEMP2=0. DO 
DO 20 1=2,10 
A1=A1+LDCI) 




B1=B1+T0(n 

TEMP1=TEriP1+TDC n^TDCl ) 

20 TEr.P2=TEnP2+TDCI)'^LD(I) 

B1=(TEMP2-TD( mAl-LD( l)=«(B1-9.0QnDn)))/ 
. (TEnpi-2.DonD( m'Bi+g.DO^TDcn^TDcn) 
A1=LDC1)-B1'*TD(1) 


A2=0.D0 


B2=0.00 


TEMP 1=0. DO 
TEMP2=0. DO 


DO 30 1=636, 6-10 
A2=A2+LD(I ) 

B2=B2+TDCI ) 

TEMP1=TEMP1+TD(n*TD(I) 

30 TEMP2=TEMP2+TDn)^L0(I) 

B2=(TEMP2-TD(64n-<A2“LD(6‘ini'(B2-5. DOnDCOmn )/ 
. CTEMP1~2. D0'*-T0(6A1 ) '*B2+5. DOl-TDCe-ll )) 
A2=LD(6-n)-B2'*-TD(6‘n) 

FIRST=, FALSE. 

GOTO 110 

200 FLUn=10. DQ-fh-CAI + BI^LOGT) 


RETURN 

300 FLUM=1Q. D0^'T:(A2+B2t^L0GT) 


RETURN 

400 WR1TE(6,6D01) 

6001 FORMATC1H ,'ODPS - WE SHOULD NOT BE HERE') 
STOP 
END 

REAL FUNCTION FG*SCS) 


THIS FUNCTION CALCULATES THE FORCE OF GRAVITY ALONG THE 
FLUX TUBE AS A FUNCTION OF S, THE DISTANCE ABOVE THE 
SOLAR SURFACE. 

IMPLICIT REALMS CA-H,0-2) 

REAL'^S RSUN/6. 9599D 1 0/, D/6 . 670-8/, MSUN/ 1 . 9S9D33/ 

LOGICAL NOTIST/. FALSE./ 

IF(N0T1ST)G0T0 10 

GM=MSUN^G 

N0T1ST=.TRUE. 

10 R=(RSUN+S) 

FG=GM/(R^-'*^2) 

RETURN 

END 

REAL FUNCTION FA*8(S) 

THIS FUNCTION CALCULATES THE AREA OF THE FLUX TUBE 
AS A FUNCTION OF S, THE DISTANCE ABOVE THE SURFACE OF 
THE SUN. 

IMPLICIT REAL'»8 (A-H.O-2) 

COMMON beta , EBETA , B 1 3, TEMP 1 , TEMPO , TCHI , D 
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REAL*8 RSUN/6. 9599D10/,A0/1 .DO/ 

LOGICAL N0T1ST/. FALSE./ 

COMMON /PARAM/ FRAC , AMP , SCALE 

IF(NOT1ST)GOTO 10 

D=FRAC*^RSUN 

R=S+D 

AR=r*»-'2 

AR=A0/AR 

FA=AR^ (R**2) 

N0T1ST=.TRUE. 

RETURN 
10 R=D+S 

FA=AR-'^(R**2) 

RETURN 

END 

REAL FUNCTION FDADS^'SCS) 

THIS FUNCTION CALCULATES THE LOGARITHMIC DERIVATIVE 
OF THE AREA AS A FUNCITON OF S, THE DISTANCE ABOVE THE 
SURFACE OF THE SUN. 

IMPLICIT REAL^^S (A-H,0-Z) 

COMMON BETA, EBETA, B 13, TEMPI ,TEMPC,TCH 1,0 
FDADS=2. D0/(D+S) 

RETURN 

END 

REAL FUNCTION FSOR^’SCS) 

THIS FUNCTION CALCULATES THE (AD HOC) NON-THERMAL 
ENERGY INPUT INTO THE SOLAR PLASMA AS A FUNCTION OF S, 

THE DISTANCE ABOVE THE SUN'S SURFACE. 

IMPLICIT REALMS (A-H,0-Z) 

REALMS RSUN/6. 9599D10/, 

.PIBY2/2-11 1921FB5AA‘120 18/ 

LOGICAL NQT1ST/. FALSE./ 

COMMON /PARAM/ FRAC , AMP, SCALE 

IF(NOT1ST)GOTO 10 

SCALE=SCALE+RSUN 

ARG=1 .DO/SCALE 

AMP=AMP'*ARG 

AMP=AMP+AMP 

ARG=ARG'^PIBY2 

SQ=S 

NOT 1ST=. TRUE. 

1 n 9R =9— 

IFCSR. GT.SCALE)GOTO 20 
C = 0C0S(ARG'<-'SR) 

FSOR=AMP'*'C*C 

RETURN 

20 FSOR=O.DO 
RETURN 
END 

SUBROUTINE AITKEN (F, X, M, ERR,XTAB, FTAB, N , SORT , EXTRAP, F INT , 
.X0IF,URK, *,*,'•«:) 

SUBROUTINE AITKEN INTERPOLATES TO FIND THE VALUE OF THE FUNCTION 
(F) AT THE POINT X. IF THE ROUTINE DOES NOT ACHIEVE THE DESIRED 
RELATIVE ERROR (ERR) USING M POINTS OR IF ROUND OFF ERROR APPEARS 
TO BE PRESENT, THE ROUTINE RETURNS THE CURRENT ERROR ESTIMATE IN 
ERR, RETURNING TO THE MAIN PROGRAM AT THE FIRST STATEMENT NUMBER 
IN THE ARGUMENT LIST. THE ROUTINE REQUIRES THE TABULATED VALUES 
IN FTAB TO BE IN ORDER OF INCREASING VALUE OF X (IN XTAB). IF 
SORT IS TRUE ON ENTRY, BOTH TABLES ARE SORTED (SEE NOTE). IF THE 
VALUE OF X IS OUTSIDE THE RANGE OF THE TABLES SUPPLIED, THE 
ROUTINE RETURNS TO THE SECOND STATEMENT NUMBER IN THE ARGUMENT 
LIST - UNLESS EXTRAP IS TRUE. IF THE ROUTINE DISCOVERS TlOO 
IDENTICAL VALUES OF X IN XTAB, THE ROUTINE RETURNS TO THE THIRD 
STATEMENT NUMBER IN THE ARGUMENT LIST. 
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It M IS GREATER THAN N OR LESS THAN 2, IT IS SET TO 10. 

IF ERR IS LESS THAN 16*’«-5, IT IS SET TO 16**-5. 

ARGUMENTS (OTHER THAN STATEMENT NUMBERS): 

F INTERPOLATED VALUE OF FUNCTION AT X (REAL - OUTPUT) 

X VALUE OF INDEPENT VARIABLE (REAL - INPUT) 

M LARGEST NUMBER OF DATA POINTS TO BE USED (INTEGER - INPUT) 

ERR REQUESTED RELATIVE ERROR (REAL - INPUT) 

XTAB TABLE OF X VALUES AT WHICH F(X) IS TABULATED 
(REAL ARRAY - INPUT) 

FTAB TABLE OF F(X) AT THE CORRESPONDING POINTS IN XTAB 
(REAL ARRAY - INPUT) 

N THE LENGTH OF TABLES XTAB AND FTAB (INTEGER - INPUT) 

SORT DETERMINES WHETHER OR NOT THE INTERNAL SORING ROUTINE 
IS TO BE USED (LOGICAL - INPUT/OUTPUT) 

EXTRAP DETERMINES WHETHER OR NOT EXTRAPOLATION OUTISDE THE 
RANGE OF THE TABLES IS ALLOWED ( LO G I CAL -INPUT) 

PINT ARRAY OF SUCESSIVE INTERPOLANTS - WORKING ARRAY 
(REAL ARRAY DIMENSION > OR = M) 

XDIF ARRAY OF DIFFERENCES BETWEEN THE POINTS AT WHICH F(X) 

IS TABULATED AND X - WORKING ARRAY (REAL ARRAY 
DIMENSION > OR = M) 

WRK WORKING ARRAY FOR CURRENT LEVEL OF INTERPOLATION 
(REAL ARRAY DIMENSION > OR = M) 

INTERNAL VARIABLES-' 

TEMP TEMPORARY STARAGE LOCATION FOR INTERMEDIATE RESULTS 

FDIFF1 PREVIOUS ABSOLUTE RELATIVE DIFFERENCE BETWEEN 

INTERPOLANTS - COMPARED WITH FDIFF2 TO CHECK FOR 
CONVERGENCE (ROUND-OFF ERROR INDICATOR) 

FDIFF2 PRESENT ABSOLUTE RELATIVE DIFFERENCE BETWEEN 

INTERPOLANTS - USED TO CHECK FOR CONVERGENCE AT 
CURRENT LEVEL (ALSO SEE FDIFF1 ABOVE) 

DIFFMAX LARGEST REPRESENTABLE FLOATING POINT NUMBER (IBM 360) 

lUP USED AS POINTER IN SORT AND INTERPOLATION 

IMID USED AS POINTER IN SORT 

IDN USED AS POINTER IN SORT AND INTERPOLATION 

XUPDIF DIFFERENCE BETWEEN X AND CLOSEST UNUSED LARGER VALUE 
IN XTAB 

XDNDIF DIFFERENCE BETWEEN X AND CLOSEST UNUSED SMALLER VALUE 
IN XTAB 

LEVEL CURRENT LEVEL OF AITKEN TRIANGULAR SCHEME 
ISTEP COUNTER FOR INTERMEDIATE INTERPOLANT LOOP 
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DONE LOGICAL FLAG TO INDICATE CURRENT LEVEL OF SHELL SORT 
IS COMPLETE 

IDISP CURRENT EXCHANGE INTERVAL IN SHELL SORT 
ILAST N MINUS IDISP - UPPER LIMIT FOR SORT DO LOOP 
I COUNTER IN SORT DO LOOP 
REMARKS: 

THE ROUTINE AS PRESENTLY WRITTEN WILL NOT WORK IN WATFIV. TO 
MAKE THE ROUTINE COMPATABLE WITH WATFIV, THREE CHANGES MUST BE 
MADE. FIRST, THE ARRAYS FINT, XDIF AND WRK SHOULD HAVE 
DIMENSION M AND THE ARRAYS XTAB AND FTAB SHOULD HAVE THE 
DIMENSION N. SECOND, THE VARIABLES M AND N SHOULD BE REMOVED 
FROM THE INTEGER DECLARATION STATEMENT. THIRD, THE STATEMENT 
WHICH CHANGES M TO 10 IF CERTAIN CONDITIONS ARE MET SHOULD BE 
DELETED. 

NOTE: 

SORT METHOD USED IS SHELL SORT - THIS METHOD MAY BE VERY 
, INEFFICIENT WHEN XTAB IS PARTIALLY SORTED. 


DECLARE VARIABLES 

REAL F ,X,ERR,XTAB( 1), FTABC 1 ) , F I NT ( 1) , XD I F ( 1 ) ,WRK( 1 ) , EPS , XUPD IF , 
.XDNDIF,FDIFF,FDIFF2, DFIMAX,TEMP 
INTEGER M,N, I STEP, ILAST, LEVEL, IDISP, I UP , I DN , I MI D , I 
LOGICAL SORT, EXTRAP, DONE 

INITIALIZE VARIABLES 

DATA EPS/Z3C 100000/, D I FMAX/Z7FFFFFFF/ 

CHECK TO SEE IF M > N OR IF M < 2, IF SO SET M TO 10 
(THIS CARO MUST BE REMOVED FOR WATFIV EXECUTION AND THE 
WORKING ARRAYS DIMENSIONED TO M) 

IF(M.LT.2..0R.M.GT.N)M=10 

CHECK TO SEE IF ERR < 16:F*-5 if SO SET IT TO 16=^*-5 
I^(ERR.LT.EPS)ERR=EPS 

CHECK TO SEE IF TABLES ARE TO BE SORTED - IF NOT GO AROUND SORT 
SECTION. 

IF(.NOT.SORT)GOTO 200 

**** SORTING SECTION BEGIN 

IDISP=N 

101 IDISP=(I0ISP+1 )/2 
ILAST=N-IDISP 

102 D0NE=.TRUE. 

DO' 103 1 = 1 , ILAST 

IFCXTABd) .LT.XTABCI+TDISP) )G0T0 103 
IFCXTABd) .EQ. XTAB (I + ID ISP) ) RETURN 3 
TEMP=.XTAB(I) 

XTABCi )=XTABd + IDISP) 

XTABd + IDISP)=TEMP 
TEMP=FTABd) 

FTABCI )=FTAB(I+IDISP) 

FTAB(I+IDISP)=TEMP 
D0NE=. FALSE. 

103 CONTINUE 
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1F(.N0T.D0NE)G0T0 102 
IFdDISP. GT. DQOTO 101 

SORTING SECTION END 

200 CONTINUE 

CHECK TO SEE IF X IS WITHIN RANGE OF TABLE - IF NOT AND IF EXTRAP 
IS FALSE RETURN TO SECOND STATEMENT IN ARGUMENT LIST 

IF(X.GE.XTAB(1))G0T0 201 

X IS BELOW LOWEST X VALUE IN XTAB - EXIT UNLESS EXTRAP IS TRUE 
IFt.NOT.EXTRAPIRETURH 2 

EXTRAP IS TRUE - SET UP POINTERS AND GO TO AITKEN 
INTERPOLATION SECTION 

IUP=2 
IDNsI 
GOTO *100 

201 CONTINUE 

CHECK TO SEE IF X IS LARGER THAN LARGEST X VALUE IN XTAB - IF NOT 
BRANCH TO SEARCH SECTION 

IF(X.LE.XTAB(N)IG0T0 300 

X IS ABOVE HIGHEST X VALUE IN XTAB - EXIT UNLESS EXTRAP IS TRUE 
1F(.N0T.EXTRAP)RETURN 2 

EXTRAP IS TRUE - SET UP POINTERS AND GO TO AITKEN 
INTERPOLATION SECTION 

lUP-N 
IONsN-1 
GOTO -lOO 

SEARCH SECTION - FIND XTAB VALUES THAT BRACKET X - USE BISECTION 

300 CONTINUE 

SET UP POINTERS FOR BISECTION 

ILIPSN 

IMID=N/2 

IDNsI 

CHECK TO SEE WHICH SIDE OF EXTABCIMIDl X IS ON AND UPDATE lUP, 
IMID AND IDN - WHEN NEW INID EOUALS IDN WE ARE DONE 

301 IFCX.GT.XTABUMIDIICOTO 302 

X LE XTABUMIDI SO IUP<=1MID 0 IMtO<s(IUP+l DN)/2 
lUP-lMID 

IMID = (IUP-HDN)/2 

IF IMID > IDN WE AREN'T DONE YET - GO BACK AND CHECK AGAIN 
OTHERWISE GO TO AITKEN INTERPOLATION SECTION 

IF(IMIDvGT,IDN)GOTO 301 - 

GOTO 'IQO 

302 CONTINUE 

X > XTAB(IMID) SO IDN<=IMID S lMID<=(IUP+lDN)/2 
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IDN=iniD 

IMID=(IUP+IDN)/2 

IF iniD > IDN WE AREN'T DONE YET - 60 BACK AND CHECK AGAIN 
OTHERWISE ENTER AITKEN INTERPOLATION SECTION 

IF(IMID.GT.IDN)G0T0 301 

END OF SEARCH SECTION 

AITKEN INTERPOLATION SECTION 

400 CONTINUE 

lUP AND IDN POINT TO FIRST TWO FUNCTION VALUES USED IN 
INTERPOLATION - INITIALIZE VARIABLES 

FDIFF2=DIFMAX 

XDNDIF=XTABtIUP)-X 

XUPDIF=XTAB(IDN)-X 

START AITKEN INTERPOLATION 

DO 401 LEVEL=1,M 

DECIDE WHICH OF THE TWO TABLE VALUES POINTED TO BY lUP AND IDN 
IS TO BE USED NEXT - THE ONE WITH XTAB CLOSER TO X 

IFCABSCXUPDIFI . GT . ABS ( XDND I F ) ) GOTO 402 

WE WILL USE lUP - PUT INFORMATION IN WORKING ARRAYS 

WRK(1 )=FTA8(IUP) 

XDIF(LEVEL)=XUPDIF 

CHECK TO SEE IF WE JUST USED THE LARGEST VALUE OF X iN XTAB 
IF SO GO TO 403 AND DO FIX UP - IF NOT UPDATE lUP AND XUPDIF 

IF(IUP.GE.N)G0T0 403 
IUP=IUP+1 

XUPDIF=XTAB(IUP)-X 

BRANCH AROUND CODE TO INTERPOLATION LOOP FOR THIS LEVEL 
GOTO 404 

FIX UP FOR USE OF LARGEST X IS TO SET XUPDIF TO LARGEST 
REPRESENTABLE FLOATING POINT NUMBER 

403 XUPDIF=DIFMAX 

BRANCH AROUND CODE TO INTERPOLATION LOOP FOR THIS LEVEL 
GOTO 404 

WE WILL USE IDN - PUT INFORMATION IN WORKING ARRAYS 

402 WRK( n=FTAB(IDN) 

XDIF(LEVEL)=XDNDIF 

CHECK TO SEE IF WE USED THE SMALLEST VALUE OF X IN X IN XTAB 
IF SO GO TO 405 AND DO FIX UP - IF NOT UPDATE IDN AND XDNDIF 

IFdDN.EQ, DGOTO 405 
IDN=IDN-1 

XDNDIF=XTAB(IDN)-X 

BRANCH AROUND CODE TO INTERPOLATION LOOP FOR THIS LEVEL 
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OOTO -10^) 


FIX UP FOR USE OF SMALLEST X IS TO SET XDNDIF TO LARGEST 
REPRESENTABLE FLOATING POINT NUMBER 

^05 X0NDIF5DIFMAX 

SKIP INTERPOLATION CALCULATION IF LEVEL IS 1 

R04 IFCLEVEL.LE, DGQTO -106 

AITKEN INTERPOLATION LOOP 

DO R07 ISTEP=2. LEVEL 
TEMP."XDIF( LEVEL) “XaiFCISTEP-n 

CHECK TO SEE IF !IE ARE GOING TO DIVIDE BY Q IF SO RETURN 
TO THIRD STATEMENT NUMBER IN ARGUMENT LIST 

I F CTEMP. EO. 0, ) RETURN 3 

CALCULATE INTERMEDIATE INTERPOLANTS 

HQ7 lIRKaSTEP)-CFINTCISTEP“1)^'XDlF(LEVEU - 

1 14RKaSTEP“l)^XDIF(lSTEP-n)/TEMP 

ENTER INTERPOLANT IN PINT 

HOB FINTCLEVEUslIRKCLEVEL) 

SKIP CHECK FDR CONVERGENCE FOR LEVEL LESS THAN H 

IFCLEVEL.LT.HIGOTO HOT 

CHECK FOR CONVERGENCE AT THIS LEVEL - IF SO BRANCH OUT 

FDirF2“i.-'»-ABS((FINT(LEVEL3-FlNTC).EVEL-1))/ 

1 (FINKLEVED + FINTCLEVEL- 1 ) ) ) 

IF(FDIFF£.LI,ERR)G0TD HOS 

SKIP ROUND OFF ERROR CHECK FOR LEVEL LESS THAN 6 
IF(LEVEL.LT.6)G0T0 HOI 

IF INTERPOLANTS ARE NOT CONVERGING - EXIT 

IF(FDIFF2,GT.FtDFFt)GOTO SOI 

UPDATE FDIFFt AND CONTINUE 
HOI FDIFF1SF0IFFS 

IF INTERPOLATED TO LEVEL-M WITHOUT CONVERGENCENCE ^ EXIT 
GOTO 501 

SET F EQUAL TO FINTCLEVEU AND RETURN 

HQ8 F=FINTCLEVEL) 

RETURN 

TERMINATIONS DUE TO LACK OF CONVERGENCE OR ROUND OFF ERROR 

sol LEVEL=LEVEL'“1 
ERR=FIDFF1 
F^iFINTCLEVEL) 

RETURN 1 
END 
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STEADY STATE ATMOSPHERE MODEL MAIN ROUTINE 


IMPLICIT REAL*8 (A-H,0-Z) 

THIS PROGRAM CALCULATES THE RUN OF TEMPERATURE, DENSITY 
HEAT FLUX AND VELOCITY IN AN INDIVIDUAL FLUX TUBE. 

THE PROGRAM READS IN PARAMETERS THAT CONTROL THE NUMBER OF 

SETS OF TABLES READ IN (NTAB), AND THE INDEPENDENT VARIABLE 

THAT CONTROLS THE FREQUENCY OF TABULATION (ITEST). FOR EACH 

SET OF TABLES THE PROGRAM READS IN THE NUMBER OF DIFFERENT INITIAL 

CONDITIONS FOR UHICH THE INTEGRATION IS TO BE PERFORMED (NRUN) 

AND A VARIABLE THAT CONTROLS WHETHER OR NOT THE RESULTS OF 
THE INTEGRATION ARE ONLY PRINTED OUT OR BOTH PRINTED OUT 
AND WRITTEN OUT IN A FORMAT SUITABLE FOR REREADING 
BY ANOTHER PROGRAM (NOUTI. IF NOUT IS LESS THAN 1, THEN THE 
RESULTS ARE ONLY PRINTED. IF NOUT IS GREATER THAN OR EQUAL TO 
1, THEN THE RESULTS OF THE INTEGRATION ARE BOTH PRINTED OUT AND 
WRITTEN OUT TO LOGICAL UNIT 10. 

THE PROGRAM READS IN A TABULATED FUNCTIONS OF DISTANCE AND 
A TABULATED FUNCTIONS OF TEMPERATURE: 

FUNCTIONS OF S= 

G THE FORCE OF GRAVITY ALONG THE TUBE 

DA THE LOGARITHMIC DERIVATIVE OF THE AREA OF THE TUBE 
WITH RESPECT TO DISTANCE ALONG THE TUBE (1/A DA/DS) 

SO A PHENOMENOLOGICAL NON-THERMAL HEAT SOURCE 

A THE AREA OF THE FLUX TUBE 

FUNCTIONS OF T: 

OK INVERSE OF THE THERMAL CONDUCTIVITY 
LU THE RADIATIVE LOSS FUNCTION 

CH THE FRACTION OF HYDROGEN NUCLEI THAT ARE IONIZED 

DC DERIVATIVE OF THE FRACITIONAL IONIZATION (CH) 

FOR EACH RUN WITH A SET OF TABLES, THE PROGRAM READS IN SO, THE 
STARTING DISTANCE, DSO THE INITIAL STEP SIZE, PRCT, THE 
MULTIPLICATIVE FACTOR BY WHICH THE INDEPENDENT VARIABLE SELECTED 
BY ITEST IS ALLOWED TO CHANGE BETWEEN TABULATION POINTS, THE 
INITIAL TEMPERATURE TO, INITIAL DENSITY NO, INITIAL HEAT FLUX QO, 
INITIAL VELOCITY UO, TSTOP, THE TEMPERATURE AT WHICH THE 
INTEGRATION WILL STOP, SSTOP, THE DISTANCE AT WHICH THE 
INTEGRATION WILL STOP, EPS, THE MAXIMUM RELATIVE ERROR IN AN 
INDEPENDENT VARIABLE ALLOWED PER DSO, AND MSTOP, THE MACH NUMBER 
AT WHICH THE INTEGRATION WILL STOP. IF THE INITIAL HEAT FLUX READ 
IN IS GREATER THAN 10 TO THE 50TH (A VERY UNPHYSICAL VALUE) THE 
INITIAL HEAT FLUX IS DETERMINED BY THE CONDITION THAT THE NET 
ENERGY FLUX IS ZERO AT THE STARTING POINT. 

THE CURRENT VERSION INTERPOLATES THE RESULTS OF THE 
INTEGRATION TO PRINT OUT VALUES OF DISTANCE TEMPERATURE, 

DENSITY, HEAT FLUX, VELOCITY, PRESSURE AND A QUANTITY WHICH 
CAN BE INFERRED FROM EUV OBSERVATIONS (P*»:2 KAPPA/Q WHERE KAPPA 
IS THE THERMAL CONDUCTIVITY) AT VALUES OF THE TEMPERATURE 
INITIALIZED IN THE ARRAY TMPOUT. IN ADDITION THE INITIAL AND 
FINAL POINTS OF THE INTEGRATION ARE PRINTED OUT. IF THE 
RESULTS ARE TO BE WRITTEN TO LOGICAL UNIT 10 ALL THE TABULATED 
RESULTS ARE PRINTED AS CALCULATED BY ABMINT. 
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DECLARE AND INITIALISE ARRAYS AND VARIABLES 

REALMS EPS,TSTOP,SSTOP,AO,UO,Nt),TO,DO,SO,S,DSO,DS,PRCT, 
.P,YPASS(‘n.KAY/ 4339 F 2 CQ 60000000 Q/* 

.MSTDP 

REAL-^S G($1D),OA(819),SO(S19),A(81S).DKCS19)iLU(S19) i CHC819) , 
.DC(819) .YTABCSi 204S) , TABU) , GRAV . DADS, SDR, AR, ORAP, LAM, CIU, DCHl 
INTEGER'^-! NRUM, NTAB , I RUN i ITAB , t , J 

EQUIVALENCE (TAB C 1) , GRAV ) , (TABC2) , DADS), (TAB ( 3) , SDR ) , (TAB(4) ,AR) , 
. CTAB(5),0KAP) , (TAB( 6) , LAM) , (TAB(7) , CHI ) , (TAB(S) ,DCHI) 

EXTERNAL DIVE 

REAL-^-S SPRIN , SB0TT/Z‘*.79 1 ‘13ECIOOOOOOOO/, TMPOUT C 'U'))/ 

, KD3, 1-503,2,03,3^3,-1.03,5.03.6.03,7.03,8^3,9.03, 

. 1 .0-1, 1.50 -1,2. OH. 3.0-1, T. 09,5.0-1, 6. DH, 7. D-1, 8. D>1, 9. 09, 

. 1.05, 1.50 5, 2. 05, 3. 05, 9, 05, 5. 05,6.05, 7.05,8.05,9.05, 

. 1 ,06, 1 . 506,2.0 6,3.06,9.06,5.06.6,06,7.06, 8.0 6,9.06/ 

5001 F0RMAT(2I5) 

5002 FORMAT (32 16) 

5003 FORMAT (90 13, 6) 

6001 FORMATdH ,T9, 'S(CM) ',T19, 'TEMP'.T39, 'N',T99, 'Q(CGS) 

• T69, 'U(CGS) ' ,T79. 'P(CGS) ' ,T99, 'LLIM' ,/) 

6002 FORMAT(7( 1PD15.6)) 

6009 F0RrnT( nil, T9. 'S(GM) ',T19, 'TEMP' ,T39, 'N',T99, 'Q(CGS) ' , 
.T69,'U(CGS) ' ,T79, 'P(CGS) ' ,T99. 'LUM',/) 

9009 FORMAT (10AS) 

CALL DIVTNT - PASS BASE ADDRESSES OF INTERPOLATION TABLES 
TO DIVF 

CALL DIVINT(G,DA,S0,A,QK,LU,GH,DC) 

READ IN NUMBER OF SETS OF TABLES AND INDEX OF INDEPENDENT 
VARIABLE THAT CONTROLS TABULATION FREQUENCY 

READC5, SDODNTAB.ITEST 
PRCT=1.05D0 

READ IN INTERPOLATION TABLES 

DO 999 ITAB=1,NTAB 
READ(5, 50D1)NRUN,NQUT 
READC9, 9009)0 
READ (9, 9009) DA 
READ (9. 9009) SO 
READ(9. 9009)A 
READ (9, 9009 )OK 
REA0(9.9QO9)LU 
REA0(9.9009)CH 
REA0(9, 9009)00 

READ IN INITIAL CONDITIONS FOR RUNS 

DO 99 IRUN=1,NRUN 
READ(S, 5002)S0.DS0,PRCT 
READ(8. 5003)T0,NO.QD.Ufl 
R E A 0 ( 8 . 50 0 3 ) T ST 0 P , SST 0 P , EPS . MSTO P 

INITIALIZE AO AND QO IF NECESSARY 

YPASS(1)=T0 

YPAS$(2)=N0 

CALL DIVFtSOjYPASS, TAB) 

AO=AR'*-NO'*-UO 

n s 1 

IF(QQ.LE.1.E50)60T0 5 

QD=-. 5DD'tU0'»-NQ'^(UQHlQ'F1 . 673S2D-29+5. D0'*KAY*T0'^( 1 .+CHI ) ) 

5 YPASS(3)=Q0 
YPASS(9)=U0 
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s=so 

DS=DS0 

NMAX=2048 


CALL ABMINT TO INTEGRATE EQUATIONS 

CALL ABMINT(S,YPASS,DIVF,DS,EPS,TSTOP,SSTOP,MSTOP, YTABC 1, II), 
.AO, ITEST, PRCT.NMAX) 
n=I 1+NMAX-1 
20 CONTINUE 

PRINT OUT RESULTS 

WE LOOK FOR VALUES OF TEMPERATURE THAT BRACKET VALUES 
OF TEMPERATURE IN TMPOUT AND INTERPOLATE. WE ALSO 
DO OUR OWN PAGINATION. 


WRITEC6, 60041 

1TEMP=0 

ILINE=1 

1=2 

SFRIN=YTAB(5, 1 ) 

CALL DIVF(SPRIN,YTAB(1, n,TAB) 

P=YTAB(2, 1 )*U . DO+CHn^YTABU , 1 ) 

RADOUT=1.D50 

IF(YTAB(3, n.EQ.O.DOlGOTO 30 
RAD0UT=-(P»-'P)/(0KAP*YTAB(3, D) 

30 P=KAY*P 

SPRIN=SPRIN-SB0TT 

WRITE! 6, 6002)SPRIN, YTABC 1, 1), YTABC 2, 1),YTABC3, 1), 
. YTABC4, 1 ) , P,RAD0UT 

FIND NEXT OUTPUT TEMPERATURE 

105 ITEMP=ITEMP+1 

IFCITEMP.GT.40)GOTO 130 

IFCYTABC 1 , I) . GT. TMPOUT C IT EMP)) GOTO 105 

FIND PRINT TEMPERATURE AND PRINT 
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115 


125 


IFCYTABC 1 , 1+1 ). GT.TMPOUTCITEMP))GOTO 115 
1 = 1 + 1 

IFC I . GE. I DGOTO 130 
GOTO 110 

FRAC=CTMP0UTCITEMP)-YTABC 1 ,I))/CYTAB( 1 , 1+1) -YTABC 1 , I) ) 
YPASSC 1)=TMP0UTCITEMP) 

YPASSC2)=YTABC2, I) + FRAC'+C YTAB C 2, 1 + D-YTABC2, I) ) 
YPASSC3)=YTABC3, I )+FRAC*C YTA3 C 3, 1+1) -YTABC 3, I ) ) 
YPASSC4)=YTABC4, D + FRAC+-CYTABC4, 1 + D-YTABC4, D) 
SPRIN=YTABC5, I ) + FRAC* C YTAB C 5 , 1 + 1 ) -YTAB C 5 , I ) ) 

CALL D I VFC YTABC 5, I), YTABC 1,1), TAB) 

P1=YTABC2, I )*C 1 . DO+CHI )*YTABC 1 , I ) 

CALL D I VFC YTABC 5, 1+1), YTABC 1, 1+1), TAB) 

P2=YTABC2, 1+1 )*C 1 . DO+CHI )^YTABC 1 , 1+1 ) 

P=P1+FRAC*CP2-P1 ) 

RA00UT=1.D50 

IFCYPASSC2) .EQ.O.DO)GOTO 125 
RAD0UT=-CP+'P)/C0KAP^'YPASSC3) ) 

P=KAY^P 

SPRIN-SPRIN^SBOTT 

WRI TEC 6, 6002) SPR IN, TMPOUT Cl TEMP) , YPASSC 2) , YPASSC 3) , 
YPASSC4) , P,RADOUT 
1=1 + 1 

IFCr.GE. I DGOTO 130 

ILINE=ILINE+1 

IFCILINE. LT. 58)GOTO 105 

ILINE=0 

WRITEC6, 6004) 

GOTO 105 
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130 IFd.GT. I DGOTO 99 
SPRIN=YTAB(5, n 

CALL DIVF(SPRIN,YTAB(1,I),TAB) 

P=YTAB(2, I )^n . DO+CKD-mABL 1, 1) 

RAD0UT=1. :D50 

IF(YTAB(3> I) .EQ.O.DO)GOTO 135 
RAD0UT=-(P'*-P)/(0KAP*YTAB(3, 1 ) ) 

135 P=KAY^P 

SPRIN=SPRIN-SBOTT 

14RITE(6,6002JSPRIN,YTAB(1,I),YTAB(2, 1)>YTAB(3, 1), 

. YTABC*!, I ) ,P, RAGOUT 
99 CONTINUE 
C ‘ 

C WRITE OUT TABULATED RESULTS OF INTEGRATION IF REQUESTED 

C 

IFCNOUT.GE. nCALL WRTR (YTAB, NMAX) 

999 CONTINUE 
STOP 
END 

SUBROUTINE WRTR(Y,N) 

IMPLICIT REALM’S (A-H,0-Z) 

REALM'S Y(5,NI 

1001 FORMATCIOAS) 

1002 FO.RMAT(A-I) 

WRITEC10, 1002IN 
WRITEC10, 100DY 
RETURN 

end _ 

DRTGINAE PA(5^ ti 
M £0QR QUAUT5C 
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ABMINT 


SUBROUTINE ABMI NT (S, Y INT , F , DS, EPS, TSTOP, SSTOP, MSTOP, YTAB, AO, 

. ITS,PRCT,NMAX) 

ABMINT SOLVES A SET OF THREE COUPLED ORDINARY DIFFERENTIAL 
EQUATIONS PLUS A CONSERVATION RELATION THAT DESCRIBE THE 
(STEADY STATE) BEHAVIOR OF A COMPRESSIBLE FLUID IN A FLUX TUBE. 

THE ROUTINE TAKE THE FOLLOMING INPUT PARAMETERS: 

S THE INITIAL DISTANCE (ARBITRARY) 

YINT THE INITIAL VALUES OF Y(1)-Y(4), THE INDEPENDENT 
VARIABLES (TEMPERATURE, DENSITY, HEAT FLUX AND 
VELOCITY) 

F THE NAME OF THE SUBROUTINE THAT CALCULATES THE 

DERIVATIVES OF THE INDEPENDENT VARIABLE AND THE VELOCITY 
(MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE 
CALLING ROUTINE) 

DS THE INITIAL STEP SIZE 

EPS THE DESIRED ACCURACY (RELATIVE) FOR A DISTANCE DS 

TSTOP THE MAXIMUM (OR MINIMUM) TEMPERATURE TO WHICH THE 
ROUTINE WILL INTEGRATE 

SSTOP THE MAXIMUM (MINIMUM) DISTANCE TO WHICH THE ROUTINE 
WILL INTEGRATE 

MSTOP THE MAXIMUM MACH NUMBER TO WHICH THE ROUTINE WILL INTEGRATE 

YTAB AN ARRAY IN WHICH THE RESULTS OF THE INTEGRATION ARE 

RETURNED TO THE CALLING PROGRAM - SHOULD BE DEMINISIONED AT ’ 
LEAST 5*NMAX. VARIABLES STORED IN THE FOLLOWING ORDER: 
TEMPERATURE, DENSITY, HEAT FLUX, VELOCI TY AND DISTANCE 

AO THE AREA AT THE STARTING POINT TIMES THE DENSITY AT THE 

STARTING POINT TIMES THE VELOCITY AT THE STARTING POINT ' 

(A CONSERVED QUANTITY) 

ITS INDEX OF THE VARIABLE THAT CONTROLS THE FREQUENCY 
AT WHICH RESULTS ARE PUT IN YTAB 

PRCT THE MULTIPLICATIVE FACTOR BY WHICH THE ITS ELEMENT OF Y 

IS ALLOWED TO CHANGE BETWEEN THE TABULATION OF THE RESULTS 

NMAX THE MAXIMUM NUMBER OF TABULATION POINTS 


THE ROUTINE USES SEVERAL LOCAL WORKING ARRAYS 
RUNGE-KUTTA: 

Y(4),Y1(4)»FO(4),Fl(‘n,F2(4) USED TO STORE INTERMEDIATE 
VALUES OF THE INDEPENDENT VARIABLE AND THEIR DERIVATIVES 

ADAMS-BASHFORTH-MOULTON PREDICTOR CORRECTOR: 

YW(32),FW(32) USED TO STORE LAST 8 STEPS OF INTEGRATION. 
THE PRESENT INTEGRATION USES 4 PREVIOUS VALUES TO ESTIMATE 
THE NEXT VALUE SO DOUBLING THE STEP SIZE CAN BE DONE IF 
AT LEAST 4 INTEGRATION STEPS HAVE OCCURRED SINCE THE LAST 
DOUBLING OF THE STEP SIZE 

YP(4),YC(4) USED TO STORE THE PREDICTED AND CORRECTED 
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i VALUES OF THE INDEPENDENT VARIABLES 

I 

AB^IINT USES AN ADAMS-BASHFORTH-MOULTON PREDICTOR-CORRECTOR 
IlNTEGRATION SCHEHE. START UP IS ACCOMPLISHED BY BACKWARD 
jlNTEGRATION WITHA RUNGE-KUTTA SCHEME AND MISSING VALUES 
NEEDED WHEN HALVING THE STEP SIZE A_RE PROVIDED USING THE 

Same runge-kutta scheme 

DECLARE VARIABLES 
IMPLICIT REALMS (A-H,0-Z) 

REAL*8 S,EPS,ERR,DS,DT,TST0P,TDIF1,TDIF2,H2,H3,H6,H8,H,T,ERR1, 
.YP(4),YC(4),YINT(4),YTAB(5,1),F0(4),F1(4),F2(4),Y1(4),Y(4),FP(4) 
.OK24/Z3FAAAAAAAAAAAAAB/,ERST,FW(32)/32*0. DO/, YW(32)/32*0.D0/, 
.MST0P,MACH,SST0P,SDIF1,SDIF2,H924,CCC1/Z4161C71C71C71C72/, 
.CCC2/Z4168E3SE3SE38E39/,CCC3/Z4141C71C71D71C72/,FRAC 
INTEGER*4 I , J,K, IND, ini , IN2, IN3, IN4, IT, DOUBLE 
L0GICAL*4 DONE 

START UP USING INTEGRATION BY 4TH ORDER R-K E 1/32 DS 


VCH1=PRCT*YINT(ITS) 

VCH2=PRCT*VCH1 

FRAC=1.D0 

iMACH=l . 649959D8*MST0P*MST0P 
MMAX=NMAX 

15 DT=.03125D0*DS*FRAC 
H2fDT*.5D0 
H3=DT/3.D0 
H6=H3*.5D0 
iH8=H2*.25D0 
T=S 
IND = 4 

DOl 111 = 1,4 
11 YTABCI, n=YINTCn 
YTAB(5,n=S 

CALL F(T,YINT,FO,AO,C999) 

DOS 1 1=1 , 4 
Y(I)=YINT(I) 

YWm=YINT(n 
1 ^ FW(I)=FO(I) 


OT Glffl C D PAGE IS 
m ?00R quality 


DO 2 1=1,3 
, DO 3 J=1,2 

' DO 4 K= 1,3 I 

4: Y1 (K)=F0(K)^H3+Y(K) 

: CALL F(T+H3,Y1,F1,A0,E999) 

DO 5 K=1,3 

5 Y1 (K)=CF0(K)+F1 (K) )*H6+Y(K) 

CALL F(T+H3,Y1,F1,A0,E999) 

DO 6 K=1,3 

6, Y1(K) = (F1(K):*:3.D0+F0(K))*H8+Y(K) 

CALL F(T+H2,Y1,F2,A0,E999) 

! DO 7 K=1,3 

7: Y1 (K)=(F2(K)*4.D0-F1(K)*3.D0+F0(K))«H2+Y(K) 

T=T+DT 

CALL F(T,Y1,F1,A0,£999) 

! DO 9 K=1 , 3 

9| Y(K) = (F2(K)»^4.D0+F1(K)+F0(K))*H6+Y(K) 

31 CALL F(T,Y,F0,A0,eg99) 

DO 8 J= 1,4 
^ YW(IND+J)=Y(J) 

8i FW(IND+J)=FO(J) 

2i IND=IND+4 

IF(YU( 12+ITS) .LT.VCHDGOTO 16 
, FRAC=FRAC*0. 500 
GOTO 15 

16 SDIF1=S-SST0P 
S=T 
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TDlF1i=YW(1)-TST0P 
H=DS*FRAC*0. 0625D0*0N24 

DT=DS»:FRAC*. 0625D0 „ 1 

H9?4=DT*.375D0 

ER|l=EPS*FRAC*.015625D0 

ERR1=ERR*0. 03125D0 

IN1=0 

IN2=4 

IN3=8 

IN4=12 

D0UBLE=4 

DONE=. FALSE. 

C 

C END INITIALIZATION -- 

C 

I=2n , I I 

50 OOUBLE=DOUBLE-1 

55 CONtiNUE 

DO 20 J=1, 3 

20 I YP( J)=Y!I(IN4+J) + H924^(CCC1*FW(IN4+J)-CCC2*FM(IN3+J) 

+ CCC3'^FU(IN2+J)-FW(IN1+J)) 

CALL F(S+DT,YP,FP,A0,£ig99) 

DO 30 J=K 3 ” 

30 i YC(J)=YW(IN4+J)+H-T(9.*FP(J) + 19.*FW(IN4+J)-5.*FW(1N3+J) 

i+FU(IN2+J)) 

TDIF2=YC( n-TSTOP 

' IF(DSIGNCTDIF2,TDIF1I ,NE.TDIF2)D0NE=.TRUE. 

SDIF2=S-SST0P 

IF(DSIGN(SDIF2,SDIF1) .NE.SDIF2)D0NE=.TRUE. 

35 CALL F(S+DT,YC,FP,A0,G999) 

IF(YC(4)^YC(4). GT.MACH*YC( 1))D0NE=.TRUE. 

ERST = OABS(YP( 1) - YC ( 1)) / ( DABS ( YP ( 1) ) + D ABS ( YC ( 1 ) ) ) 

; ERST=DMAX1 (DABS(YP(2)-YC(2))/(DABS(YP(2) )+DABS(YC(2))) ,ERST) 

I IF(DABS(YC(3) )+0ABS(YP(3) ) . LT. 1 . E-30IGOT0 220 

ERST = DMAXUDABS(YP(3)-YC(3) J/(DABS(YP(3) ) + DABS(YC(3.))),ERSTJ 
220 IF(0ABS(YC(4))+DABS(YPC4)).LT. 1.E-30)GOTO 215 

ERST-DMAX1 (DABS(YP(4)-YC(4) )/(0ABS(YP(4))+DABS(YC(4)) ) ,ERST) 
215 IF(YC(ITSKGT.VCH2)G0T0 90 

IFCERST.GT.ERRIGOTO 100 
IFIERST . LT. ERRDGOTO 200 
205 S=S+DT 

IN1=IN2 
IN2=IN3 
-TN3=IN4 

IN4=M0DCIN4+4,32) 

DO 40 J= 1,4 

YW(IN4+J)=YC(J) 

40 FU(IN4+J)=FP(J) 

IF(YC(ITS).LT.VCH1)G0T0 50 
DO 60 J=1 , 4 

60 ; YTABC J, I)=YC(J) 

YTAB(5, I)=S 
1 = 1+1 

IFCI.GT.NMAXIGOTO 10 
VCH1=VCH2 
; VCH2 = PRCT-TVCH1 

I IF(DONE)GOTO 10 . ^- - - 

! GOTO 50 , 

90 IACC=1 

GOTO 103 
100 IACC=0 
’ 103 CONTINUE 
C i 

C SECTION THAT HALVES INTEGRATION STEP 
C ^ 

IFCDT.LT. 1 .D-DGOTO 205 
IFrOT.GT.DSlGOTO 109 
ERR=ERR*. 5D0 
ERR1=ERR1*.5D0 
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109 DT=DT*.5D0 
H=DT»:0N24 
H924=DT*.375D0 

REARRANGE WORKING ARRAYS 

IT=M0D(IN4+12,32) 

J=1 N'l 

DO' i01 K=1,4 

FU(IT+K)=FW(IN4+K) 

101 YW(IT+K)=YW(IN4+K) 
TN4=IT 

IT=MOD(IN3+8, 32) 

DO 102 K=1,4 

FW(IT+K)=FW(IN3+K) 

YW(IT+K)=YW(IN3+K) 

FW(IN3+K)=FU(IN2+!<) 

102 VW(IN3+K)=YU(IN2+K) 
IN2=IT 

IT=IN3 

IN1=J 

lN3=M0D(IN2+4,32) 


Pj6 JOOB. 





generate missing information with 4TH ORDER R-K 


H2=DT*.25D0 

H3=DT/6.D0 

'' 'R6-=H3*.5D0 

H8=H2*. 2500 
T=S-(DT+DT) 

DO 110 J= 1,4 
! Y(J)=YW(IN2+J) 

lilO F0(J)=FW(IN2+J) 

DO 120 J=1 ,2 
! DO 130 K=1,3 

1|30 Y1 (K)=F0(K)*H3+Y(K) 

I CALL F(T+H3,Y1,F1,A0,£999) 

00 140 K= 1,3 I 

1'40 Y1(K) = (F0(K) + F1(K))*H6+Y(K) 

CALL F(T+H3,Y1,F1,A0,C999) 

DO 150 K=1,3 

150 Y1 (K) = (F1 (K)*3.D0+F0(K))*H8+Y(K) 

CALL F(T+H2,Y1,F2,A0,£999) 

DO 160 K=1 , 3 

160 YKK) = (F2(K)*4.D0-F1(K)*3.D0+F0(K))*H2+Y(K) 

; T=T+H2+H2 

CALL F(T,Y1,F1,A0,S999) 

'T- D0180K=1,3 

180 Y(K)=(F2(K)*4.00+F1(K)+F0(K))*H6+Y(K) 

CALL F(T, Y,F0,A0,£999) 

120 1 CONTINUE 

DO 170 J= 1,4 

FU(IN3+J)=F0(J) ' ; 

170 YW(IN3+J)=Y(J) ' ’ 

: DO 115 J= 1,4 - 

Y(J)=YW(IT+J) 

115 F0(J)=FW(IT+J) - 

T=S-4. D0»-DT 
i DO 125 J=1,2 
! DO 1135 K=1,3 

135 Y1(K)=F0(K)*H3+Y(K) 

CALL F(T+H3,Y1,F1,A0,£999) 

DO 145 K=1, 3 

145 YKK) = (F0(K) + F1(K))*H6+Y(K) 

CALL F(T+H3,Y1,F1,A0,£999) 

DO 155 K=1,3 

155 Y1(K)=(F1(K)*3.D0+F0(K))*H8+Y(K) 

CALL F(T+H2,Y1,F2,A0,£999) 

DO 165 K=1,3 
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'-"JiV 

165 Y1(K)=(F2CK)*4.D0-F1(K)*3.D0+F0(K))*H2+Y(K) 

T=T+H2+H2 

CALL F(T,Y1,F1,A0,C999) 
i DO 185 K=1,3 

185 Y(K)=(F2CK)*4.D0+FT(X)+F0(K))*H6+Y(K) 

! CALL F(T, Y, F0,A0,£999) 

125 CONTINUE 

DO 175 J=1,‘1 

FWCIN1+J)=F0(J) 

175 YU(IN1+J)=YCJ) 

GOTO 55 

RETURN TO ABM P-C INTEGRATION WITH NEW STEP SIZE 
200 CONTINUE 

SECTION THAT DOUBLES INTEGRATION STEP SIZE 

If (!ACC. EQ. 1 )GOTO 205 
IF (DOUBLE. GE.OIGOTO 205 
00;UBLE=‘1 
S^S-DT 
.:‘UT = 0T+DT 
■'H = DT'^0N24 

H92*1 = DT*.375D0 . ' 

IF(0T.GT.DS)G0T0 209 

ERR=ERR+ERR 

ERRI-ERRI+ERRI 

209 K=M0BCrN4+4,32) 

IT=M0D(:iN4+12, 32) 

DO 210 J=1 

FWilN^I+JI^FWdNB+J) 

YW(INT+J)=YW(IN3+J) 

FW(IN3+J)=FW(IN1+J) I 
YW(IN3+J)=YW(IN1+J) : 

FW(IN1+J)=FW(K+J) 
i YW( IN1+J)=YW(K+J) 

I FW(IN2+J)=FW(IT+J) 

210 YW( 1N2+J)=YW( IT+J) 

GOTO 205 

10 NMA>! = I-1 
RETURN 

999 CONTINUE 

6001 FORMATCIH, 'FATAL ERROR T OR S OUT OF TABULATED RANGE') 
WRITE(6,6001) 

RETURN 
END , 
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DRIUINXC TO5E B 
OP POOR QUALITY 


OIVF 


1 

I 


I CSECT 


01VjF(S,Y,DY,A0,*) OR D I VF ( S, Y , TAB) 

REALMS S,Y(4),DY(4),A0 REALMS S, Y ( 4) , TAB ( 8 ) 

i 

THd FIRST FORM OF THE CALL CALCUALTES THE DERIVATIVES 

OT/DS = DY(1)/DS , DN/DS = DY(2)/DS t DO/DS = DY(3)/DS 

AND STORES THEN IN ARRAY DY. THE VALUE OF V=YC4) IS COMPUTED 
FROM THE CONSERVATION LAW NVA = CONSTANT AND STORED IN Y(4). 
THE ROUTINE INTERPOLATES THE VALUES OF PRARMETERS NEEDED FOR 
THE CALCULATIONS FROM TABULATIONS OF 4 FUNCTIONS OF S ONLY AND 
4 FUNCTIONS OF T ONLY. IF S OR T IS OUT OF THE TABULATED 
RANGE, THE OFFENDING QUANTITY IS NEGATED AND THE ROUTINE DOES 
the equivalent of a FORTRAN RETURN1. 

THE SECOND FORM OF THE CALL STATEMENT (DISTINGUISHED FROM THE 
FIRST BY THE NUMBER OF ARGUMENTS) CALCULATES THE INTERPOLATED 
! values of g(S),da/ds(s),source(s), area(s) and I/KAPPACT), 
■laHbdact), CHICT) and dchi/dtct) and stores them in tab. 

;THERE is a second ENTRY POINT (DIVINT) WHICH PICKS UP 
;AN0 STORES LOCALLY THE ADDRESSES OF THE TABULATIONS OF 
|THE FUNCTIONS NEEDED FOR THE CALCULATIONS. 

'NOTE THAT THIS MEANS THAT MEANINGLESS RESUI-TS 
WILL BE PRODUCED IF DIVINT IS NOT CALLED BEFORE 
THE FIRST TIME DIVF IS CALLED. IT IS EVEN POSSIBLE 
THAT SOME SORT OF ABEND WILL RESULT. 


ITHE following two FORTRAN SUBROUTINES ARE ROUGHLY EQUIVALENT 
TO THE TWO CALLS TO DIVF (DIVINT IS NOT REPRODUCED) 

i SUBROUTINE D I VF (S, Y , DY , AO, *) 

IMPLICIT REALMS (A-H.O-Z) 

DIMENSION Y(4),DY(4) 

REALT8 KOMH/S. 29S9776D7/,TION/I . 0464606D 5/, KAY/1 .380620-16/ 
DY(1)=-Y(3)*QKAP(T) 

Y(4)-A0/(Y(2)»'A(S) ) 

DY(2) = (Y(2)*(Y(4)*Y(4)*DA0S(S)-DY( 1 )=*^KOMH 
. ( ( 1 . DO + CHI (T) ) + (T'*^0CHI (T) )-G(S)*S) ) )/ 

. (KOMH’^d . DO+CHI (T))‘«Yn)-Y(4)»Y(4) ) 

DY (3)=-(L(T)*CHI (T)^Y(2)*Y(2))-( 1 . 5D0*Y (4)*Y ( 2)*’KAY*DY ( 1 ) 

> .'^■((I .DO+CHI (T)) + (Y(1)-TION)»-DCHL(T)))+((1 .DO_CHI (T))'*KAY 

.*Y(1)^Y(4)^DY(2))+S0R(S)-Y(3)*DADS 
RETURN 
: END 

i 

i SUBROUTINE 0 1 VF (S , Y , TAB) 

! IMPLICIT REALMS (A-H,0-Z) 

: DIMENSION Y(4),TAB(8) 

TAB(T) G(S) ! 

TABL2)=DADS(S) ; 1 
TABC3)=S0R(S) i 

TABC4)=A(S) ' 

TAB(5)=OKAP(T) 

TAB(6)=L(T) 

TABC7) = CHI (T) . 

TABC8)=DCHI (T) 

RETURN 

END 

NOTE: IN THE COMMENTS 'R' RE:FERS TO GENERAL PURPOSE REGISTERS 
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IVINT 


FIRST 


DFIRST 


* 


* 


* 


* 

* 

* 

* 

* 



AND 'F' REFERS TO FLOATING POINT REGISTERS. 

USING *,15 TELL ASSEMBLER NEST INST ADDR IN R15 

B DFIRST BRANCH AROUND NAME AND OTHER ENTRY POINT 

DC X'04‘ 

DC CLB'DIVF ' 

DIVINT(G,DADS,SOR,AREA,OKAP,LUM,CHI,DCHI) 

REAL*8 G ( 563), DADS ( 563), SOR( 563), AREA (563), 

REAL*8 0KAP(820),LUM(820),CHI(820) ,DCHI (820) 

ENTRY DIVINT 


USING 

*,15 

TELL ASSEMBLER NEXT INSTR ADDR IN R15 

B 

TFIRST 

BRANCH AROUND NAME 

DC 

X '06 ' 


DC 

CL7'DIVINT ' 


STM 

14, 12, 12(13) 

SAVE CALLING ROUTINES GPR'S 

LM 

2, 9, 0(1) 

GET BASE ADDR'S OF INTERPOLATION TABLES 

STM 

2,9,GADDR 

SAV TABLE BASE ADDR'S 

LM 

2,9,28(13) 

RESTORE CALLING ROUTINE'S GPR'S 

MV I 

12(13), X'FF' 

INDICATE CONTROL RETURNED 

BR 

14 

RETURN FROM INITIALIZATION 


MAIN ROUTINE RESUMES 
USING DIVF, 1 5 

STM H, 12, 12(13) SAVE CALLING ROUTINES GPR'S 
L 2,0(1) R2 <= ADDR S 

LM 3,6,GADDR R3-R6 <= ADDR'S OF TABLES FOR S 

MVC FL0AT+1(6),2(2) FLOAT <= FRACTIONAL DISTANCE FROM 

NEXT SMALLER VALUE OF S TABULATED. 


LH 

12,0(2) 

R12 <= HIGH ORDER BYTES OF S 

S 

12,SDISP 

REDUCE R12 BY SDISP - « WORDS FROM 
BASE OF TABLES 

LD 

4, FLOAT 

F4 <= FRACTION 0 LE FRAC LE 1 

BM 

BADS 

IF RESULT IS NEGATIVE OUT OF RANGE 
GOTO BADS 

C 

12,SBND 

COMPARE R12 TO SBND - IF GREATER 

BH 

BADS 

OUT OF RANGE GOTO BADS 

SD 

4,=D' .5' 

F4 <= X = FRAC - .5 -.5 LE X LE .5 

SLA 

12,3 

R12 <= R12*8 NOW BYTE DISPLACEMENT 
FROM BASE OF INTERPOLATION TABLE. 

NOW COMPUTE WEIGHTS 
FUNCTIONS OF S 

FOR CUBIC INTERPOALTION OF 

LDR 

2,4 

F2 <= X 

MDR 

4,4 

F4 <= X**2 = X2 

HDR 

4,4 

F4 <= X2/2 

SD 

4,=D' 1. 125' 

F4 <= X2/2 - 9/8 

LDR 

6,4 

F6 <= X2/2 - 9/8 

HDR 

4,4 

F4 <= X2/4 - 9/16 

MDR 

6,2 

F6 <= X3/2 - 9X/8 

LCDR 

0,4 

FO <= -X2/4 + 9/16 

ADR 

0,6 

FO <= X3/2 - X2/4 - 9X/8 + 9/16 

STD 

0,WM1 

WM1 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 
ING TO CLOSEST SMALLER VALUE OF S. 

LCDR 

0,4 

FO <= -X2/4 + 9/16 

SDR 

0,6 

FO <= -X3/2 - X2/4 + 9X/8 + 9/16 

STD 

0,WP1 

WP1 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 
ING TO CLOSEST LARGER VALUE OF S. 

ADR 

6,2 

F6 <= X3/2 - X/8 

AD 

4,=D'.5' 

F4 <= X2/4 - 1/16 

MD 

6, =X' 4055555555555555' F6 <= X3/6 - X/24 

LDR 

0,4 

FO <= X2/4 - 1/16 

SDR 

0, 6 

FO <= -X3/6 + X2/4 + X/24 - 1/16 

ADR 

6,4 

F6 <= X3/6 + X2/4 - X/24 - 1/16 

STD 

0,WM3 

WM3 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 
ING TO 2ND CLOSEST SMALLER VALUE OF S. 
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STD 6,WP3 


WP3 <= WEIGHT FOR TABLE ENTRY CORRESPOND- 
ING TO 2ND CLOSEST LARGER VALUE OF S. 


NOW CALCULATE INTERPOLATED VALUES OF GRAVITY AND DA/DS 
(HAVE WEIGHTS FOR TABLES ENTRIES 1 S 4 IN FPR'S 0 C 6) 


LDR 

LDR 

MD 

MO 

MO 

MD 

ADR 

ADR 

LD 

LDR 

MD 

MD 

ADR 

ADR 

LO 

LDR 

MD 

MO 

ADR 

ADR 

STD 

STD 

LD 

LD 


4, 
2 , 
0 , 
2 , 
4, 
6 , 
0 , 
4, 
2 , 
6 > 
2 , 
6 , 
0 , 
4, 
2 , 
6» 
2 . 
6 , 
0 , 
4, 
0 , 
4. 
0 , 
2 , 


0 

6 

0(3. 12) 
24(3, 12) 
0(4,12) 
24(4,12) 
2 
6 

WM1 

2 

8(3,12) 
8(4, 12) 

2 

6 

WP1 

2 

16(3, 12) 
16(4, 12) 
2 
6 
G 

DADS 

WN3 

WP3 


F4 < = 
F2 < = 
FO < = 
F2 < = 
F4 < = 
F6 < = 
FO < = 
F4 < = 
F2 < = 
F6 < = 
F2 < = 
F6 < = 
FO < = 
F4 < = 
F2 < = 
F6 < = 
F2 < 
F6 < 
F2 < = 
F4 < = 
G < = 
DAOS 
FO < = 
F2 < = 


GRAV 
GRAV ' 
DA/DS 
DA/DS 


PSGE B 
OP. POOR QUALITY 


WEIGHT 
WEIGHT 
WEIGHT 
WEIGHT 
WEIGHT 
WEIGHT 

W1*G1 + W4*G4 
Wr^DI + W4^D4 
WEIGHT 2 
WEIGHT 2 
W2*G2 
W2*D2 

W1*G1 + W4-"^G4 
W1*D1 + W4*D4 
WEIGHT 3 
WEIGHT 3 
= U3*G3 
= W3=*-D3 

INTERPOLATED VALUE OF G 
INTERPOLATED VALUE OF DA/OS 
INTERPOLATED VALUE OF GRAVITY 
<= INTERPOLATED VALUE OF DA/DS 
WEIGHT 1 
WEIGHT 4 


W2*G2 

W2*D2 


* NOW CALCULATE 

* ( HAVE WEIGHTS 

* 


VALUES OF SOURCE AND AREA 
1 8 4 IN FPR'S 0 £ 2) 


LDR 

4,0 

F4 

<= WEIGHT 1 

LDR 

6,2 

F6 

<= WEIGHT 4 

MD 

0,0(5, 12) 

FO 

<= WEIGHT 1 * SOURCE 1 

MD 

2,24(5, 12) 

F2 

<= WEIGHT 4 * SOURCE 4 

MD 

4,0(6,12) 

F4 

<= WEIGHT 1 * AREA 1 

MD 

6,24(6,12) 

F6 

<= WEIGHT 4 area 4 

ADR 

0,2 

FO 

<= W1*S1 + W4*S4 

ADR 

4,6 

F4 

<= W1=»-A1 + W4»^A4 

LD 

2,WM1 

F2 

<= WEIGHT 2 

LDR 

6,2 

F6 

<= WEIGHT 2 

MD 

2,8(5, 12) 

F2 

<= W2‘«S2 

MD 

6, 8(6, 12) 

F6 

<= W2S-A2 

ADR 

0.2 

FO 

<= W1*S1 + W4*S4 + W2"FS2 

ADR 

4,6 

F4 

<= W1*A1 + W4»-'A4 + W2-'*-A2 

LD 

2.WP1 

F2 

<= WEIGHT 3 

LDR 

6,2 

F6 

<= WEIGHT 3 

MD 

2, 16(5, 12) 

F2 

: <= W3^S3 

MD 

6, 16(6, 12) 

F6 

. <= W3^-'A3 

ADR 

0,2 

F2 

<= INTERPOLATED VALUE OF 

ADR 

4,6 

F4 

<= INTERPOLATED VALUE OF 

STD 

0,S0R 

SOR 

<= INTERPOLATED VALUE Oi 

STD 

4, AREA 

AREA <= INTERPOLATED VALUE I 


SOURCE 
AREA 

F SOURCE 
OF AREA 

* 

* CALCULATE INDEX AND FRACTIONAL DISPLACEMENT FOR INTERPOLATION 

* ON TEMPERATRUE (T) TABLE 


2,4(1) 


LM 

3, 6,KA0DR 

MVC 

FL0AT+1 (6 

LH 

12,0(2) 

S 

12,TDISP 

LD 

4, FLOAT 


R2 <= BASE ADDR Y ARRAY 
R3-R6 <= BASE ADDR'S TABLES FOR 
INTERPOLATION OF FUNCTIONS OF T 
U FLOAT <= FRACTIONAL DISTANCE FROM 
NEXT SMALLER VALUE OF T TABULATED. 
R12 <= HIGH ORDER BVTES OF T 
REDUCE R12 BY TDISP - # WORDS FROM 
BASE OF TABLES 

F4 <= FRACTION 0 LE FRAC LE 1 
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BM 

BADS 

IF RESULT IS NEGATIVE OUT OF RANGE 
GOTO BADS 

C 

12,TBND 

COMPARE R12 TO TBND - IF GREATER 

BH 

BADS 

OUT OF RANGE GOTO BADS 

SD 

4,=D' .5' 

F4 <= X = FRAC - .5 -.5 LE X LE .5 

SLA 

12,3 

R12 <= R12*8 NOW BYTE DISPLACEMENT 
FROM BASE OF INTERPOLATION TABLE. 

NOW COMPUTE WEIGHTS 
FUNCTIONS OF T 

FOR CUBIC INTERPOALTION OF 

LDR 

2,4 

F2 <= X 

MDR 

4,4 

F4 <= X**^ = X2 

HOR 

4,4 

F4 <= X2/2 

SD 

, 4,=D' 1. 125' 

F4 <= X2/2 - 9/8 

LDR 

6,4 

F6 <= X2/2 - 9/8 

HDR 

4,4 

F4 <= X2/4 - 9/16 

MDR 

6,2 

F6 <= X3/2 - 9X/8 

LCDR 

0,4 

FO <= -X2/4 + 9/16 

ADR 

0,6 

FO <= X3/2 - X2/4 - 9X/8 + 9/16 

STD 

0,WM1 

UM1 <= WEIGHT FOR TABLE ENTRY CORRESPOND 
ING TO CLOSEST SMALLER VALUE OF T. 

LCDR 

0,4 

FO <= -X2/4 + 9/16 

SDR 

0,6 

FO <= -X3/2 - X2/4 + 9X/8 + 9/16 

STD 

0,WP1 

WP1 <= WEIGHT FOR TABLE ENTRY CORRESPOND 
ING TO CLOSEST LARGER VALUE OF T. 

ADR 

6,2 

F6 <- X3/2 - X/8 

AD 

4,=D' ,5' 

F4 <= X2/4 - 1/16 

MO 

6, =K '4055555555555555' F6 <= K3/6 - X/24 

LDR 

0,4 

FO <= X2/4 - 1/16 

SDR 

0,6 

FO <= -X3/6 + X2/4 + X/24 - 1/16 

ADR 

6 , 4 

F6 <= X3/6 + X2/4 - X/24 - 1/16 

STD 

0,WM3 

WM3 <= WEIGHT FOR TABLE ENTRY CORRESPOND 
ING TO 2ND CLOSEST SMALLER VALUE OF T. 

STD 

6,WP3 

UP3 <= WEIGHT FOR TABLE ENTRY CORRESPOND 
ING TO 2ND CLOSEST LARGER VALUE OF T. 

NOW CALCULATE INTERPOLATED 

VALUES OF 1/KAPPA(T) AND L 

(HAVE WEIGHTS FOR TABLES ENTRIES 1 E IN FPR'S 0 E 6) 

LDR 

4,0 

F4 <= WEIGHT 1 

LDR 

2,6 

F2 <= WEIGHT 4 

MD 

0,0(3, 12) 

FO <= WEIGHT 1 * K 1 

MD 

2,24(3,12) 

F2 <= WEIGHT 4 * K 4 

MO 

4,0(4, 12) 

F4 <= WEIGHT 1 * L 1 

MD 

6,24(4,12) 

F6 <= WEIGHT 4 * L 4 

ADR 

0,2 

FO <= W1*K1 + W4^FK4 

ADR 

4,6 

F4 <= W1*L1 + W4*L4 

LD 

2,WM1 

F2 <= WEIGHT 2 

LDR 

6,2 

F6 <= WEIGHT 2 

MD 

2, 8(3,- 12) 

F2 <= WZ'i-KZ 

MO 

6,8(4,12) 

F6 <= W2^L2 

ADR 

0,2 

FO <= W1*K1 + W4^K4 + W2*K2 

ADR 

4,6 

F4 <= Wl^ll + W4*L4 + W2'«L2 

LD 

2,WP1 

F2 <= WEIGHT 3 

LDR 

6,2 

F6 <= WEIGHT 3 

MD 

2,16(3,12) 

F2 <= WS'FKS 

MD 

6,16(4,12) 

F6 <= W3*L3 

ADR 

0,2 

F2 <= INTERPOLATED VALUE OF 1/KAPPA 

ADR 

4,6 

F4 <= INTERPOLATED VALUE OF L 

STD 

0,0KAP 

OKAP <= INTERPOLATED VALUE OF 1/KAPPA 

STD 

•4, LUM 

LUM <= INTERPOLATED VALUE OF L 


* CHECK TO SEE IF GAS FULLY IONIZED (T > 65,536, BYTE OISP > 38C(HEX)), 
IF SO SKIP INTERPOLATION OF CHI AND DCHI/OT, IF NOT CONTINUE 

C 12, MBYTE 

BC 2, MIGHT 

LO 0,WM3 FO <= WEIGHT 1 
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LO 


2.WP3 


F2 <= WEIGHT 4 


* NOW CALCULATE VALUES OF CHI AND DCHI/DT 

* ( HAVE WEIGHTS 1 £ 4 IN FPR'S 0 £ 2) 


LDR 

4,0 

F4 <= WEIGHT 1 


LDR 

6,2 

F6 <= WEIGHT 4 


MO 

0,0(5, 12) 

FO <= WEIGHT 1 * CHI 1 


MD 

2,24(5, 12) 

F2 <= WEIGHT 4 * CHI 4 


MD 

4,0(6, 12) 

F4 <= WEIGHT 1 * DCHI/DT 1 


MD 

6,24(6, 12) 

F6 <= WEIGHT 4 * DCHI/DT 4 


ADR 

0,2 

FO <= W1*C1 + W4*C1 


ADR 

4,6 

F4 <= WIi^DI + W4*D4 


LD 

2,WM1 

F2 <= WEIGHT 2 


LOR 

6,2 

F6 <= WEIGHT 2 


MD 

2,8(5, 12) 

F2 <= W2*C2 


MD 

6,8(6, 12) 

F6 <= W2^'02 


ADR 

0,2 

FO <= W1*C1 + W4*C4 + W2^C2 


ADR 

4,6 

F4 <= W1=^D1 + W4*D4 + W2^:D2 


LD 

2,WP1 

F2 <= WEIGHT 3 


LDR 

6,2 

F6 <= WEIGHT 3 


MD 

2, 16(5, 12) 

F2 <= W3*C3 


MD 

6, 16(6, 12) 

F6 <= W3*D3 


ADR 

0,2 

F2 <= INTERPOLATED VALUE OF 

CHI 

ADR 

4,6 

F4 <= INTERPOLATED VALUE OF 

DCHI. 

STD 

0,CHI 

CHI <= INTERPOLATED VALUE OF 

CHI 

STD 

4,DCHI 

DCHI <= INTERPOLATED VALUE OF DC 

LM 

3, 4, 8(1) 

GPR'S 3£4 (= BASE ADDR'S OF 

DY ; 

LTR 

3,3 

CHECK IF 3 NEGATIVE IF SO IS 

LAST 

BM 

TABL 

PARAMETER - GOTO TABL 


LD 

2,0KAP 

F2 <= 1/KAPPA 


LNDR 

2,2 

F2 <= -1/KAPPA 


MD 

2,16(2,0) 

F2 <= -Q/KAPPA = DT/DS 


STD 

2, 0(3,0) 

STORE DT/DS 


AD 

0,=D' 1. ' 

FO <= 1+CHI 


LD 

2, 0(4,0) 

F2 (= AO = NO^^UO^'AO 


LD 

6, AREA 

F6 <= AREA 


MD 

6, 8(2,0) 

F6 <- AREAi^N 


DDR 

2,6 

F2 <= (NO*UO*AO)/(AREA*N) = 

U 

STD 

2,24(2,0) 

STORE U (VELOCITY) 


MDR 

2,2 

F2 <= U**2 


LDR 

6,0 

F6 <= 1+CHI 


MD 

6,K0MH 

F6 <= (K*( 1+CHI ) )/MH 


MD 

6, 0(0, 2) 

F6 <= (KT*( 1+CHI ) )/MH 


SDR 

6,2 

F6 <= (KT*( 1+CHI) )/MH - U**2 


MD • 

2, DADS 

F2 <= DADS+’U**2 


SD 

2,G 

F2 <= DADS:*-U**2 - GS 


MD 

4, 0(2,0) 

F4 <= DCHI/DT*T 


ADR 

4,0 

F4 <= ( 1+CHI ) + DCHI/DT*T 


MD 

4,K0MH 

F4 <= K/MH«F4(^) 


MD 

4, 0(3,0) 

F4 <= DT/DS»!F4('t) 


SDR 

2,4 

F2 <= F2(4^) - F4(*) 


MD 

2, 8(2,0) 

F2 <= F2(t)=^'N: 


|*(U**2*'0A/DS - DT/DS*K/MH*((1+CHn+T*DCHI/DT) - GS) 

DDR 

2,6 

F2 <= DN/DS 


STD 

2,8(3, 0) 

STORE DN/DS 


MDR 

2,0 

F2 <- ( 1+CHI)*DN/DS 


LD 

4, 0(2,0) 

F4 <= T 


MDR 

2,4 

F2 <= T*( 1+CHI )*DN/DS 


SD 

4,TI0N 

F4 <= T - TI 


MD 

4, DCHI 

F4 <= (T-TI)^'DCHI/DT 


ADR 

4,0 

F4 <= (1+CHI) + (T-TI)*DCHI/DT 

LO 

0, 8(2,0) 

FO <= N 


MDR 

4,0 

F4 <= N*( 1+CHI +(T-TI)*DCHI/DT) 

MD 

4,=D' 1.5' 

F4 <= 3/2 * F4(t) 


MD 

4, 0(3,0) 

F4 <= DT/DS * F4(*) : 






3/2*DT/DS*( 1+CHI + 

(T-TI )*DCHI/DT) 


SDR 

2,4 

F2 <= (1+CHI)*DN/DS*T - 

F4(«t) 

MD 

2,24(2,0) 

F2 <= F2(6) * U 


MD 

2, KAY 

F2 <= K * F2('t) 


NOR 

0,0 

FO <= N**2 


MD 

0, CHI 

FO <= N**2 * CHI 


MD 

0,LUM 

FO C= RADIATIVE LOSSES 


SDR 

2,0 

F2 <= F2(t) - FOIil:) 


LD 

0, 16(2,0) 

FO <= Q 


MD 

0, DADS 

FO <= Q/A*DA/DS 


SDR 

2,0 

F2 <= DQ/DS LESS SOURCE 

TERM 

AD 

2,S0R 

F2 <= DQ/DS 


STD 

2, 16(3,0) 

STORE DQ/DS 


LM 

14,12, 12(13) 

GPR'S RETURNED TO ORIGINAL STATE 

MV I 

12(13),X'FF' 

TELL CALLING PROGRAM WE 

'RE RETURN 

BCR 

15,14 

RETURN 



END OF SECTION FOR T , 65,536 NEXT SECTION DOES SAME CALCULATIONS 
FOR FULLY IONIZED CASE 


LM 

3, 4, 8(1) 

FPR'S 3S4 <= BASE ADDR'S OF DY E AO 

MVC 

0KAP+16(16),=X 

'41100000000000000000000000000000' 
CHI <= 1 £ OCHI <= 0 

LTR 

3, 3 

CHECK IF 3 NEGATIVE IF SO IS LAST 

BM 

TABL 

PARAMETER - GOTO TABL 

MD 

0, 16(2,0) 

FO <= Q/KAPPA 

LCDR 

0,0 

FO <= -Q/KAPPA = DT/DS 

STD 

0, 0(3,0) 

STORE DT/DS 

LD 

2,0(4, 0) 

F2 <= AO = NO*UO=«AO 

LD 

4, AREA 

F4 <= AREA 

MD 

4, 8(2,0) 

F4 <= AREA^N 

DDR 

2,4 

F2 <= (NO'^UO*AO)/(AREA'^-N) = U 

STD 

2,24(2,0) 

STORE U (VELOCITY) 

MDR 

2,2 

F2 <= U'*'*2 

LDR 

6,0 

F6 <= DT/DS 

LD 

4,K0MH2 

F4 <= 2K/MH 

MDR 

6,4 

F6 <= 2K/MH*DT/DS 

MD 

4, 0(2,0) 

F4 <= 2KT/MH 

SDR 

4, 2 

F4 <= 2KT/MH - U**2 

MD 

2, DADS 

FS <= U'»'''*-2/A*DA/DS 

SDR 

2,6 

F2 <= Ut=*2/A*DA/DS - 2K/MH*DT/DS 

SD 

2, G 

F2 <= F2(i) - G 

LD 

6, 8(2,0) 

F6 <= N 

MDR 

2,6 

F2 <= F2($) * N 

DDR 

2,4 

F2 <= DN/DS 

STD 

2, 8(3,0) 

STORE ON/DS 

MDR 

0,6 

FO <=N:»^DT/DS 

LDR 

4,0 

THIS AND THE NEXT TWO INSTRUCTIONS 

ADR 

0, 0 

EFFECTIVELY MULTIPLY FO BY 3. 

ADR 

0,4 

FO <= 3N*DT/DS 

MD 

2, 0(2,0) 

F2 <= DN/DS*T 

ADR 

2,2 

F2 <= 2T*DN/DS 

SDR 

2,0 

F2 <= 2T»’DN/DS - 3N*DT/DS 

MD 

2, KAY 

F2 <= K*(2T'^DN/DS - 3N*DT/DS) 

MD 

2,24(2,0) 

F2 <= K'^U'^(2T*DN/DS - 3N*DT/DS) 

MDR 

6, 6 

F6 <= N^^*2 

MD 

6, LUM 

F6 <= RADIATIVE LOSSES 

SDR 

2,6 

F2 <= F2(4) - F6(4') 

LD 

4, 16(2,0) 

F4 <= Q 

MD 

4, DADS 

F4 <= Q/A*DA/DS 

SDR 

2,4 

F2 <= DQ/DS LESS SOURCE TERM 

AD 

2,S0R 

F2 <= DQ/DS 

STD 

2, 16(3, 0) 

STORE DQ/DS 

LM 

14, 12,12(13) 

GPR'S RETURNED TO ORIGINAL STATE 

MV I 

12(13), X'FF' 

TELL CALLING PROGRAM WE'RE RETURNING 

BR 

14 

RETURN 

RETURN 

INTERPOLATED 

FUNCTION VALUES - NOT DERIVITIVES 
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« 


TABL 

me 

0(64,3),G 

PUT INTERPOLATED VALUES IN TAB 


LM 

14,12, 12(13) 

RESTORE GPR'S 


MV I 

12(13), X'FF' 

INDICATE CONTROL RETURNED 


BR 

14 

RETURN 

* CASE 

OF S OR 

T OUTSIDE OF 

THE TABULATED RANGE 

BADS 

LO 

0, 0(2,0) 

FO <= OFFENDING QUANTITY (T OR S) 


LNDR 

0,0 

FO NOM NEGATIVE 


STD 

0, 0(2,0) 

STORE OFFENDING QUANTITY - FLAG 


LM 

14, 12, 12(13) 

GPR'S RETURNED TO ORIGINAL STATE 


LA 

15,4 

GPR 14 <= 4 (RETURN 1 ) 


MV I 

12(13), X'FF' 

TELL CALLING PROGRAM ME'RE RETURNING 


BR 

14 



CNOP 

4,8 


* STORAGE FOR 

ADDR'S, CONSTANTS, AND INTERNAL VARIABLES 

GADDR 

DC 

X'OOOOOOOO' 


DADDR 

DC 

X'OOOOOOOO' 


SADDR 

DC 

X'OOOOOOOO' 


AADDR 

DC 

X'OOOOOOOO' 


KADDR 

DC 

X'OOOOOOOO' 


LADDR 

DC 

X'OOOOOOOO' 


CADDR 

DC 

X'OOOOOOOO' 


GCDDR 

DC 

X'OOOOOOOO' 


SOISP 

TDISP 

DC 

DC 

X'00004710' 

X'00004410' 

OlUGlKAL PAdfi 18 

TEND 

SBND 

DC 

DC 

X'00000330' 

X'00000330' 

QB £QQR QUAUT1C 

HBYTE 

DC 

X'00000778' 



MM3 DC X'OOOOOOOOOOOOOOOO' 

MM1 DC X'OOOOOOOOOOOOOOOO' 

MP1 DC X'OOOOOOOOOOOOOOOO' 

MP3 DC X'OOOOOOOOOOOOOOOO' 

G DC X'OOOOOOOOOOOOOOOO' 

DADS DC X'OOOOOOOOOOOOOOOO' 

SOR DC X'OOOOOOOOOOOOOOOO' 

AREA DC X'OOOOOOOOOOOOOOOO' 

OKAP DC X'OOOOOOOOOOOOOOOO' 

LUM DC X'OOOOOOOOOOOOOOOO' 

CHI DC X'OOOOOOOOOOOOOOOO' 

DCHI DC X'OOOOOOOOOOOOOOOO' 

TION DC X'45198C6100000000' 

KAY DC X'339F2CB600000000' 

KOMH DC X'474EAB1BOOOOOOOO' 

KOMH2 DC X'479D5A3600000000' 

FLOAT DC X'4000000000000000' 

END 
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